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PREFACE 


The  purpose  of  this  thesis  was  three-fold.  The  first  purpose  was 
to  revise  and  extend  the  capabilities  of  existing  software  for  selecting 
the  significant  control  variables  of  a  simulation  model,  based  on  a 
newly  developed  selection  criterion.  The  second  purpose  was  to  compare 
the  results  obtained  using  the  revised  software  employing  two  different 
selection  procedures.  And  the  third  purpose  was  then  to  validate  the 
effectiveness  of  the  new  selection  criterion  by  comparison  to  results 
derived  using  other  common  selection  criteria. 

After  extensive  revision,  the  software,  now  renamed  the  Variable 
Subset  Selection  Program  (VSSP) ,  was  ready  for  use.  The  VSSP  was  then 
used  to  evaluate  data  with  known  characteristics  and  data  derived  from 
an  untested  simulation  model.  The  results  obtained  from  this  effort 
served  to  demonstrate  the  usefulness  of  the  VSSP  and  the  validity  of  the 
the  new  selection  criterion.  It  is  highly  recommended  that  the  work  be 
continued,  as  further  benefits  are  yet  to  be  realized  and  may  be  of 
substantial  significance. 

The  execution  and  preparation  of  t  .  _  thesis  would  not  have  been 
possible  without  the  help  of  others.  I  a.  eeply  indebted  to  my  faculty 
advisor,  Major  Kenneth  Bauer,  Jr.,  for  his  extensive  time,  patience,  and 
assistance.  I  also  wish  to  thank  my  thesis  reader,  Lt  Colonel  Thomas 
Schuppe,  for  pointing  out  my  numerous  writing  errors  and  ensuring  the 
final  product  was  understandable.  Finally,  I  wish  to  thank  my  family 
and  friends  for  their  continuous  support  and  encouragement  when  the 
going  got  rough. 

James  A.  Gigliotti 
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ABSTRACT 


./  The  purpose  of  this  thesis  was  three-fold.  The  first  purpose  was 
to  revise  and  extend  the  capabilities  of  existing  software  for  selecting 
the  significant  control  variables  of  a  simulation  model,  based  on  a 
newly  developed  selection  criterion.  The  second  purpose  was  to  compare 
the  results  obtained  using  the  revised  software  employing  two  different 
selection  procedures.  And  the  third  purpose  was  then  to  validate  the 
effectiveness  of  the  new  selection  criterion  by  comparison  to  results 
derived  using  other  common  selection  criteria. 

Extensive  revision  of  the  existing  software  was  necessary  to 
prepare  it  for  use.  Initially,  the  software  was  revised  to  extend  its 
adaptability  to  evaluating  new  data  and  to  increase  user  friendliness. 
Next,  a  new  procedure  was  added  to  the  software  to  permit  it  to  evaluate 
data  using  a  Stepwise  (Forward  Selection)  procedure.  Previously,  the 
software  only  performed  analysis  of  the  data  through  an  Enumerated 
Subsets  approach.  After  revision  of  the  software  was  complete,  it  was 
renamed  the  Variable  Subset  Selection  Program  (VSSP) . 

Once  the  VSSP  was  ready,  it  was  used  to  evaluate  two  types  of  data. 
The  first  type  of  data  was  created  using  a  known  stochastic  structure. 
Three  sets  of  this  data  was  used,  each  set  using  a  different  covariance 
structure  between  the  responses  and  control  variables.  The  second  type 
of  data  was  created  from  an  untested  simulation  model  .'"'This  data  ... 
provided  a  means  of  validating  the  program  and  the  selection  criterion 
incorporated  into  it.  In  addition,  the  data  derived  from  the  untested 
simulation  model  was  also  evaluated  using  a  commercially  available 
statistical  software  package  employing  several  common  selection 


criterions.  These  results  were  then  compared  to  those  obtained  using 
the  VSSP. 

Overall,  this  study  found  that  as  the  amount  of  data  evaluated  by 
the  VSSP  increased,  any  differences  between  the  control  variables 
selected  as  significant,  by  either  the  Enumerated  Subsets  or  Stepwise 
procedure,  disappeared.  In  fact,  a  point  was  apparently  reached  where 
additional  data  caused  no  change  in  the  results  obtained.  Also,  when 
the  covariances  between  the  control  variables  are  known,  this  only  makes 
any  difference  when  a  minimal  amount  of  data  is  available.  And  finally, 
comparison  of  the  results  obtained  by  the  VSSP  and  the  commercial 
software  package  showed  the  new  criterion  to  be  comparable  to  those 
commonly  in  use.  The  new  criterion  also  had  the  advantage  of  not 
requiring  a  subjectively  determined  stopping  criteria  for  selecting  the 
significant  control  variables,  unlike  some  of  the  other  criterion  in  use 
today . 

The  recommendations  made  from  this  study  involved  further  work  on 
the  VSSP  and  additional  experimentation  which  can  be  performed  to  extend 
the  usefulness  of  the  new  criterion.  Several  suggestions  for 
enhancements  to  the  VSSP,  primarily  in  regards  to  adding  additional 
evaluation  procedures  and  increasing  program  efficiency,  were  noted. 
There  is  also  much  work  remaining  in  regards  to  the  new  selection 
criterion.  One  possibility  mentioned,  would  be  to  see  if  the  point 
where  further  data  provides  no  additional  benefit  to  the  evaluation, 
could  be  analytically  determined. 


COMPARISON  OF  SELECTION  PROCEDURES 
AND  VALIDATION  OF  CRITERION  USED  IN 
SELECTION  OF  SIGNIFICANT  CONTROL  VARIATES 
OF  A  SIMULATION  MODEL 


Background 

When  dealing  with  computer  simulations  it  is  typically  desirable  to 
have  a  general  understanding  of  how  the  simulation  inputs  will  affect 
the  final  results.  It  is  also  desirable  to  be  able  to  accurately 
estimate  the  expected  simulation  response.  Furthermore,  if  the 
estimation  of  the  response  can  be  achieved  with  a  subset  of  the 
simulation  inputs  (variables) ,  a  variance  reduction  on  the  estimator  of 
the  mean  can  also  be  realized.  One  way  of  achieving  these  goals  is 
through  the  identification  of  a  good  subset  of  control  variates. 

Control  variates,  also  known  as  control  variables,  are  variables  which 
have  a  significant  covariance  with  the  response  of  interest. 

The  development  of  a  quick  and  easy  method  for  identifying  the 

subset  of  significant  control  variates  in  a  simulation  model  would 

greatly  decrease  the  time  and  effort  required  to  gain  insights  into  the 

simulation.  Identifying  the  significant  control  variates  for  a 

simulation  model  can  also  enhance  the  process  of  preparing  and 

implementing  an  experimental  design.  It  would  eliminate  the  guesswork 

in  determining  which  variables  to  concentrate  on  in  a  subsequent 

experimental  design.  This  could  also  save  computer  time  by  identifying 

a  subset  of  the  available  control  variates  to  work  with,  since  the 

k 

standard  experimental  design  requires  2  simulation  runs  to  acquire 
data,  where  k  is  the  number  of  variables  being  tested. 
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The  need  for  research  into  this  problem  and  several  methods  for 
approaching  it  have  been  identified  in  current  literature,  but  little 
substantial  work  has  yet  been  accomplished  (Bauer,  1987:2). 

Furthermore,  Pritsker  (1986:748)  notes  that  even  though  theoretical 
development  of  control  variates  has  proceeded,  little  practical 
application  has  been  reported. 

Specific  Problem 

The  purpose  of  this  thesis  was  to  compare  selection  procedures  for 
selecting  the  significant  control  variates  of  a  simulation  model  and  to 
validate  the  selection  criterion  used. 

The  scope  of  this  thesis  was  confined  to  revising  and  adding  new 
procedures  to  previously  written  software  for  identifying  the 
significant  control  variates  of  a  simulation  model  and  applying  the 
software  to  selected  data  sets.  The  data  was  also  evaluated  using 
commercial  software  and  additional  selection  criteria.  The  results  were 
then  used  to  compare  tb<;  selection  procedures  employed  and  to  validate 
the  selection  criterion. 

Sub-objectives 

In  order  to  solve  the  specific  problem  the  following  sub-objectives 
or  steps  were  accomplished.  The  first  sub-objective  was  to  revise  and 
incorporate  new  procedures  into  existing  software  for  evaluating 
simulation  output  and  identifying  the  significant  control  variates. 

The  second  sub-objective  was  to  test  the  revised  software  on 
several  sets  of  simulation  model  output  with  known  responses, 
significant  control  variates,  and  covariances  between  the  control 
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variates.  This  output  is  referred  to  as  the  control  data/output  in 
later  text. 

The  third  sub-objective  was  to  compare  the  sets  of  significant 
control  variates  identified  by  the  revised  software  with  each  other  and 
with  those  known  for  the  control  data/output. 

The  fourth  sub-objective  was  to  use  the  revised  software  to 
identify  the  significant  control  variates  of  an  untested  simulation 
model.  This  simulation  model  is  referred  to  as  the  test  model  or  data 
generation  model  in  later  text. 

And  the  fifth  sub-objective  was  to  compare  the  selected  control 
variates  with  variable  subsets  selected  using  commercially  available 
software  and  various  other  selection  criteria. 

General  Methodology 

For  each  of  the  sub-objectives  outlined  above,  there  were 
associated  sets  of  methods  and  equipment  required  to  accomplish  them. 
Accomplishment  of  the  first  sub-objective,  to  revise  and  incorporate  new 
procedures  into  the  existing  software  for  evaluating  simulation  output 
and  identifying  the  significant  control  variates,  necessitated  a  two¬ 
fold  approach.  The  first  approach  involved  revising  existing  software 
previously  developed  by  Bauer  (1989B).  The  basis  for  the  software  was 
an  evaluation  of  all  possible  combinations  (enumerated  subsets)  of  the 
control  variates  involved.  To  successfully  perform  this  task  required  a 
thorough  understanding  of  the  underlying  logic  and  statistical  concepts 
on  which  the  software  was  based.  The  next  step  was  the  actual  revision 
of  the  existing  software  with  the  goals  of  increasing  generality  and 
user-friendliness  of  the  software. 
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The  second  approach  to  meeting  this  sub-objective  dealt  with  adding 
a  new  routine  to  the  revised  software  based  on  a  stepwise  evaluation 
procedure  and  the  same  statistical  accept/reject  criterion  as  the 
enumerated  subset  routine.  Again,  study  was  necessary  to  understand  the 
stepwise  procedure  and  construct  the  logic  for  implementation.  When 
that  step  was  completed,  the  actual  stepwise  procedure  software  was 
written,  debugged,  and  incorporated  into  the  overall  revised  program. 

The  resulting  software  product  is  referred  to  as  the  Variable  Subset 
Selection  Program  (VSSP)  in  later  text. 

The  second  sub-objective,  test  the  Variable  Subset  Selection 
Program  on  several  sets  of  simulation  model  output  with  known  responses, 
significant  control  variates,  and  covariances  between  the  control 
variates,  was  completed  as  follows.  The  first  step  was  to  obtain 
data/output  with  these  characteristics.  Next,  this  data/output  was 
evaluated  using  the  Variable  Subset  Selection  Program.  The  data/output 
was  evaluated  using  both  the  enumerated  subsets  and  stepwise  procedures 
incorporated  into  the  program.  The  end  result  of  this  sub-objective  was 
a  single  subset  of  control  variables  derived  using  each  of  the  selection 
procedures . 

The  third  sub-objective,  compare  the  sets  of  significant  control 
variates  identified  by  the  Variable  Subset  Selection  Program  with  each 
other  and  with  those  known  for  the  control  data/output,  was 
straightf orward  and  involved  answering  a  series  of  questions.  Were 
there  any  differences  in  the  number  of  control  variates  identified  as 
significant?  Were  there  any  differences  in  the  specific  control 
variates  identified  as  significant? 
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The  fourth  sub-objective,  use  the  Variable  Subset  Selection  Program 
to  identify  the  significant  control  variates  of  an  untested  simulation 
model,  was  completed  as  follows.  The  first  step  was  to  obtain  a 
simulation  model  for  evaluation.  Next,  the  simulation  model  was  run  on 
the  AFIT  VMS  computer  system  to  create  the  output  data  to  use  as  input 
data  for  the  Variable  Subset  Selection  Program.  And  finally,  the 
Variable  Subset  Selection  Program  evaluated  the  model  output  and  a 
subset  of  the  significant  control  variates  was  selected.  The  difference 
between  the  data/output  derived  from  this  model  and  the  data/output  used 
in  sub-objective  two  is  that  the  significant  control  variates  had  not 
been  previously  determined. 

The  fifth  and  final  sub-objective,  compare  the  selected  control 

variates  with  variable  subsets  selected  using  commercially  available 

software  and  various  other  selection  criteria,  was  accomplished  as 

follows.  First,  the  SAS  statistical  package,  installed  on  the  AFIT  VMS 

system,  was  selected  as  representative  of  commercially  available 

software.  Next,  the  SAS  procedures  for  Enumerated  Subsets,  and  Stepwise 

o 

(using  Forward  Selection,  Backward  Selection,  and  R  Maximization  [MAXR] 
options)  evaluation  was  applied  to  the  data/output  of  the  untested 
simulation  model.  And  finally,  the  results  obtained  using  SAS  were 
compared  to  those  obtained  using  the  Variable  Subset  Selection  Program. 
The  primary  purpose  of  this  sub-objective  was  to  validate  the  criterion 
used  by  the  VSSP  and  demonstrate  it  will  provide  comparable  results. 

Thesis  Organization  and  Development 

A  review  of  literature  relevant  to  this  thesis  is  presented  in 
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Chapter  II.  The  literature  review  covers  the  topic  of  control  variates 
and  their  theoretical  development.  Also  reviewed  are  the  common 
selection  criteria  and  selection  procedures  in  use  today,  and  the 
theoretical  development  of  a  new  selection  criterion. 

In  Chapter  III  the  detailed  methodology  used  in  approaching  and 
completing  this  thesis  is  covered.  Then,  the  results  of  the  research 
are  presented  and  discussed  in  Chapter  IV.  And  finally,  the  conclusions 
and  recommendations  reached,  after  evaluating  the  data,  are  given  in 
Chapter  V. 
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II.  LITERATURE  REVIEW 

The  following  discussion  is  a  review  of  the  literature  that  has 
relevance  to  this  thesis  topic. 

Variance  Reduction  Techniques 

The  need  for  some  form  of  Variance  Reduction  Technique  becomes 
apparent  when  it  is  understood  that  simulation  is  an  experimental 
technique,  for  analyzing  systems  which  usually  involve  the  use  of 
stochastic  processes  (Tomick , 1988 : 1 .7) .  Since  stochastic  processes  are 
"a  collection  of  random  variables'  (Ross,  1985:72),  then  the  output  from 
a  simulation  experiment  is  also  a  random  variable.  Thus,  the  response 
of  interest  is  only  an  estimate  of  the  true  value.  From  Pritsker 
(1986:742),  'the  variance  of  the  sample  mean  is  a  derived  measure  of  the 
reliability  that  can  be  predicted  if  a  simulation  experiment  is 
repeatedly  performed'.  Pritsker  also  states  that  'Variance  Reduction 
Techniques  (VRTs)  are  methods  that  attempt  to  reduce  the  estimated 
values  of  variance  through  the  setting  of  special  conditions  or  through 
the  use  of  prior  information.' 

In  a  survey  of  Variance  Reduction  Techniques  (VRTs) ,  performed  by 
Wilson  (1984:280),  VRTs  are  divided  into  two  categories:  correlation 
methods,  and  importance  methods.  His  paper  discusses  three  correlation 
methods  (common  random  numbers,  antithetic  variates,  and  control 
variates)  and  four  importance  methods  (importance  sampling,  conditional 
Monte  Carlo,  stratified  sampling,  and  systematic  sampling). 

The  basic  difference  between  the  two  categories  is  the  underlying 
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principle  of  the  methods.  The  correlation  methods  increase  the 
efficiency  of  the  simulation  by  exploiting  linear  correlations  among  the 
simulation  responses  and  input  variables.  The  importance  methods 
achieve  variance  reduction  by  concentrating  on  prior  knowledge  of  the 
input  domain. 

Of  the  VRT  methods  discussed  by  Wilson,  this  thesis  concentrates  on 
the  use  of  control  variates.  The  rationale  for  this  decision  were  two¬ 
fold.  First,  this  is  a  promising  technique  which  can  provide  valuable 
insights  into  the  problem,  if  even  to  identify  a  lack  of  correlation 
between  the  inputs  (i.e.  control  variates)  and  responses.  And  second, 
even  though  theoretical  development  has  been  proceeding,  not  much  has 
been  done  in  the  way  of  practical  applications.  Further  information  on 
this  method  follows. 

Control  Variates 

The  method  of  control  variates,  also  known  as  control  variables,  is 
one  of  the  correlation  methods  mentioned  previously.  Basically,  ‘ . .  the 
control  variates  technique  uses  regression  methods  to  exploit  any 
inherent  correlation  between  an  output  and  a  selected  random  variable 
vector  with  known  mean  that  is  observed  on  each  run’  (Wilson,  1984:280). 

The  remainder  of  the  discussion  on  control  variates  will  cover  the 
types  of  control  variates,  the  theory  behind  the  concept  of  control 
variates,  and  a  summary  of  recent  work  accomplished  on  this  topic. 

Types  of  Control  Variates.  Law  and  Kelton  (1982:359)  define  two 
types  of  control  variates.  The  two  types  are  internal,  or  concomitant, 
control  variates,  and  external  control  variates.  The  first  type  of 
control  variate  addressed  is  the  internal  control  variate.  Internal 
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control  variates  nay  be  selected  from  among  any  of  the  input  random 
variables,  or  simple  functions  of  them,  since  their  means  are  known. 

This  view  is  further  endorsed  by  Bauer  ( 1989A:0-63) ,  when  he  states  'any 
input  random  variable  is  a  candidate  for  a  control  variate.'  In 
addition,  analysis  of  the  simulation's  use  of  the  input  random  variables 
should  identify  at  least  the  sign  of  the  correlation  with  the  output 
random  variable.  An  advantage  of  internal  control  variates  is  they  must 
typically  be  generated  anyway  during  a  simulation  and  therefore  add 
essentially  nothing  to  the  cost  of  running  the  simulation. 

The  second  type  of  control  variates,  external  control  variates, 
require  the  simultaneous  simulation  of  a  similar,  but  analytically 
tractable,  system  using  common  random  numbers.  The  corresponding  output 
from  this  similar  simulation  is  then  used  as  the  control  variate.  By 
analytically  tractable,  it  is  meant  that  the  expected  value  of  the 
output  variable  can  be  calculated  exactly.  It  is  then  hoped  that  the 
similar  nature  of  the  tractable  simulation  will  induce  a  correlation 
between  the  two  outputs,  which  can  then  be  exploited.  The  major 
disadvantage  of  this  type  of  control  variate  is  that  it  requires  a 
second  simulation  model  and  additional  simulation  runs,  so  it  is  not 
costless . 

Theory  of  Control  Variates.  To  restate,  ‘the  concept  associated 
with  control  variates  is  the  identification  of  a  variable,  say  X,  that 
has  a  positive  covariance  with  the  variable  of  interest,  say  T* 
(Pritsker,  1986:748). 

Unless  noted  otherwise,  the  following  theory  is  based  on  a  class 
handout  provided  by  Bauer  (1989A). 
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Univariate  Simulation  Response  with  a  Single  Control .  Assume 


that  Y  is  an  unbiased  estimator  of  the  response  of  interest  0;  that  is, 
E(Y)  =  0,  where  E(Y)  is  the  expected  value  of  Y.  Let  X  be  an  input 
random  variable,  selected  as  the  control  variate.  It  is  further  assumed 
that  X  has  a  known  expected  value  of  ux  and  is  highly  correlated  with  Y. 
Then,  for  any  constant  b  (known  as  the  control  coefficient),  the 
controlled  estimator  Y(b),  given  by  Eq  (1),  is  unbiased  for  0. 

Y (b)  =  Y  -  b  (X  -  ux)  (1) 

Then  the  variance  of  Y(b)  is 

Var [ Y (b) 3  =  Var(Y)  +  b2Var(X)  -  2bCov(Y,X)  (2) 

From  review  of  Eq  (2) ,  it  is  readily  apparent  that  a  variance  reduction 
can  be  achieved  if 

2bCov (Y,X)  >  b2Var (X)  (3) 

So,  if  the  condition  of  Eq  (3)  is  met,  then  the  controlled  estimator 
will  have  a  smaller  variance  then  the  uncontrolled  estimator.  It  is 
also  apparent  that  if  the  variables  X  and  Y  are  independent,  in  which 
case  Cov(Y.X)  =  0,  then  no  improvement  over  the  uncontrolled  estimator 
is  possible.  Next,  with  the  application  of  some  calculus  to  Eq  (2),  the 
optimal  control  coefficient,  fi,  for  which  the  variance  of  Y(b)  is  a 
minimum,  is  given  by 


fi  =  Cov (Y,X) /Var (X)  (4) 

Substituting  Eq  (4)  into  Eq  (1)  leads  to  Eq  (5)  which  gives  the  optimal 
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controlled  estimator  Y(8) 


Y (8)  =  Y  -  tCov(Y,X)/Var(X) ]  *  (X  -  uy)  (5) 

And  substituting  Eq  (4)  into  Eq  (2)  yields  the  corresponding  minimum 
variance  for  Y(fi)  of 

Var [Y(fl) ]  =  (1  -  pxy2)  »  Var(Y)  (6) 

where  pXy  is  the  correlation  coefficient  between  X  and  Y.  Therefore  as 
the  absolute  value  of  pXy  tends  to  its  maximum  value  of  one,  the 
variance  of  Y(fi)  decreases. 

For  the  next  step,  let  6  be  denoted  by  uy.  Then  the  average  of  the 
controlled  observations  Yj,  for  i  =  1  to  K,  is  an  unbiased  point 
estimator  of  Uy.  This  estimator  is  represented  by  Eq  (7). 

K 

Y(fi)  =  ( 1/K)  £  Y^fi)  (7) 

1  =  1 

where  K  is  the  sample  size  and 

Yi(fi)  =  Yi  -  fi(Xi  -  ux)  (8) 


Typically,  the  optimal  value  &  is  unknown  and  must  be  estimated. 

However,  fi  can  be  estimated  as  follows: 

An  intuitive  estimate  of  8  replaces  the  right-hand  side  of  Eq  (4) 
with  the  appropriate  sample  quantities.  This  solution  turns  out  to 
be  the  least  squares  solution  for  8.  When  the  assumption  of  joint 
normality  between  Y  and  X  is  made,  then  the  least  squares  solution 
is  also  the  maximum  likelihood  solution.  (Bauer,  1987:6) 

So  the  following  equation  provides  an  estimate  of  8. 

K  -  -  K  -  2 
fi  =  £  (Y±  -  Y)*(Xi  -  X)/  £  (X*  -  xr  (9) 

i  =  l  i  =  1 


1  1 


where 


K 

Y  =  Z  Yi/K  (10) 

i  =  l 


and 


K 

X  =  L  X^K  (11) 

i  =  1 


The  point  estimate  of  uy  is 

Y (fi)  =  Y  -  fi(X  -  ux)  (12) 

Then,  the  variance  of  the  point  estimator  is  given  by 

VarCY(fi)]  =  Var[Y(fl)]/K  (13) 

where 

Var [ Y (fl) ]  =  (l-pxy2)*Var(Y)  (14) 

Bauer  (1087:6)  provides  the  derivation  of  the  interval  estimate 
through  the  application  of  regression  theory  and  assuming  that  Y  and  X 
are  jointly  normal  random  variables.  The  resulting  100(l-a)Z  confidence 
interval  is  given  by 

-  *  -  *  I/O 

Y(B)  +  tK.2(l-a/2)*(Var[Y(fi)  ]»sn)  (15) 

where 

K  K 

sn  =  L  (Xt  -  Ux!2/Kl  (Xi  -  X)2  (16) 

i  =  l  i*  1 

tx_2  is  the  Student’s  t-distribution  with  K-2  degrees  of  freedom,  and 
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‘a‘  is  the  desired  significance  level. 


Since  £  is  estimated,  the  variance  reduction  achieved  is  smaller 
then  could  have  been  obtained  had  the  optimal  control  coefficient  been 
known.  This  loss  of  variance  reduction  is  quantified  as  the  Loss  Factor 
(LF) .  The  loss  factor  is  defined  as  ‘the  ratio  of  the  variance  of  the 
estimator  of  Uy  when  the  optimal  control  coefficient  is  not  known  to  the 
variance  of  the  estimator  when  the  coefficient  is  known"  (Bauer, 

1987:9).  Bauer  (1987:10)  provides  the  derivation  of  the  loss  factor, 
which  reduces  to 


LF  =  (K-2) / (K-Q-2) 


(17) 


where 

Q  =  the  number  of  controls 
K  =  the  number  of  independent  replications 

Furthermore,  the  ‘loss  factor  acts  as  a  multiplier  to  the  minimum 

variance  ratio  (MVR) *  (Bauer,  1987:10),  which  is  given  by 

MVR  =  Var [Y(fi) ]/Var(Y)  (18) 

The  MVR  represents  the  possible  variance  reduction  when  the  optimal 
control  coefficient  is  known.  Multiplying  Eq  (17)  and  Eq  (18)  together 
leads  to  the  variance  ratio  (VR) .  The  VR  represents  the  possible 
variance  reduction  when  £  is  estimated. 

VR  =  LF  *  MVR  (19) 

Dnivariate  Simulation  Response  with  Multiple  Controls. 

Kleijnen  (1974:151)  addresses  the  extension  of  theory  to  multiple 
control  variates.  Also,  Bauer  provides  a  summary  of  ‘the  development 
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presented  by  Lavenberg  and  Welch  (1981)  for  simulation  output  analysis 
based  on  independent  replications,  batch  means,  and  regenerative 
analysis'  (Bauer,  1987:  11). 

2 

Let  Y  be  the  univariate  response  with  variance  vy  ,  X  be  the  Qxl 
vector  of  controls,  vXy  be  the  Qxl  vector  of  covariances  between  Y  and 
X,  and  Ex  be  the  QxQ  covariance  matrix  of  the  controls.  Then  rewriting 
Eq  (12)  with  multiple  controls  leads  to 


Y(fi)  =  Y  -  flT(X  -  ux) 


where  C,  X,  and  ux  are  Qxl  vectors.  The  vector  of  optimal  control 
coefficients  is  given  by 


fi  =  Ex  vxy 


Since  the  covariance  matrices  are  usually  unknown,  £  can  be  estimated  by 

w  A 

substituting  the  sample  values  of  Ex  and  vxy  into  Eq  (21) .  This  leads 
to  the  following  equation: 


fi  =  Sx  Sxy 


where  Sx  is  the  inverse  of  the  QxQ  sample  covariance  matrix  of  the 
controls,  and  Sxy  is  the  Qxl  vector  of  sample  covariances  between  the 
univariate  response  and  the  vector  of  controls. 

Assuming  that  Y  and  X  have  a  joint  multivariate  normal 
distribution,  then 


*  2  ** 

'y  vyx 


vxy  Ex 
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Consequently,  Y(fi)  is  unbiased  for  uy  and  an  exact  100(l-a)X  confidence 
interval  is  given  by 

Y(fl)  *  tK-Q-i (l-a/2)D»SyX  (24) 

where 

D2  =  K~ 1  ♦  (K- 1 ) " 1 (X  -  ux)TSx_1(X  -  ux)  (25) 

Syx2  =  (K-Q-l)'1 (K-l) (Sy2  -  SxyTSx_1S_1Sxy)  (26) 

tj(_Q_i  is  the  Student’s  t-distribution  with  (K-Q-l)  degrees  of  freedom, 

2 

and  Sy  is  the  sample  variance  of  Y  (Bauer  et  al ,  1988:3). 

Multiple  Simulation  Responses  with  Multiple  Controls.  Bauer, 
Venkatraman,  and  Wilson  (1987:334)  provide  the  necessary  theoretical 
structure  for  handling  the  case  of  P  response  variables  and  Q  control 
variates.  When  dealing  with  multiple  variables,  Y  is  a  Pxl  vector  of 
response  variables,  fi  is  a  PxQ  matrix  of  control  coefficients,  and  S  is 
the  sample  covariance  matrix  of  the  response  vector.  Assuming  that  Y 
and  X  have  a  joint  multivariate  normal  distribution,  then 


Y 

Uy 

Ey  Eyx 

~ 

Nq+i 

«v 

» 

~ 

X 

1 

«x 

! 

Exy  Ex 

(27) 


Consequently,  Y(fi)  is  an  unbiased  estimator  of  uy  and  an  exact  100(1- 
a/2)X  confidence  ellipsoid  for  uy  is  given  by 


[!(fi-uy]TSyx‘1[Y(fl)-uy]  <  P (K-Q-l) (K-P-Q) ~*DF(l-a;P,K-P-Q)  (28) 


where 

D2  =  K'1  ♦  (K- 1 ) ~ 1 (X  -  ux)TSx_1(X  -  ux) 


(29) 
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Syx  =  (K“Q~1)  (K"l) (Sy  -  SyxSx  Sxy)  (30) 

and  F(l-a;n^,n2)  is  the  F-distribution  with  nj  and  n2  degrees  of  freedom 
(Bauer  et  al ,  1987:335). 

'The  advantage  of  the  above  approach  over  selecting  separate 
controls  for  each  response  is  the  capability  to  form  a  joint  confidence 
region  for  the  response  vector,  rather  than  being  limited  to  univariate 
confidence  intervals'  (Tomick,  1988:2.10). 

Selection  of  Significant  Control  Variates 

Neter  states  that  ‘One  of  the  most  difficult  problems  in  regression 
analysis  often  is  the  selection  of  the  set  of  independent  variables  to 
be  employed  in  the  model*  (1983:417).  Regardless  of  the  problem 
involved,  there  are  several  reasons  to  restrict  the  number  of  variables 
used  in  a  model:  (1)  A  model  with  a  large  number  of  variables  can  be 
expensive  to  maintain,  (2)  Models  with  a  limited  number  of  variables  are 
easier  to  analyze  and  understand,  and  (3)  The  presence  of  many  highly 
intercorrelated  variables  may  add  little  to  the  predictive  power  of  the 
model,  detracting  from  the  model’s  descriptive  abilities  and  increasing 
the  problem  of  roundoff  error  (Neter,  et  al;  1983:418). 

The  selection  of  the  significant  control  variates  depends  primarily 
on  the  selection  criteria  and  selection  procedure  used.  The  selection 
criteria  determines  the  relative  significance  of  a  regression  variable 
(control  variate)  and  this  may  vary  as  the  criteria  varies.  The  type  of 
selection  procedure  has  an  effect  on  whether  the  subset (s)  of  control 
variates  chosen  is  the  ‘best*  subset  or  is  a  'near-best'  subset. 

Selection  Criteria.  The  most  common  selection  criteria  in  use 
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In  addition,  a 


2  2 

today  and  which  are  reviewed  here,  are  Rp ,  Ra ,  and  Cp. 
new  selection  criteria,  BCp,  is  also  presented.  Any  of  these  selection 
criteria  can  be  used  for  selecting  one  or  more  variable  subsets. 

In  the  discussion  of  each  criteria,  the  following  notation  is  used. 
P  is  the  number  of  potential  parameters.  The  intercept  term,  ,  counts 

as  one  parameter,  so  there  are  P-1  potential  X  variables  (X ^ . Xp_ j) . 

p  is  the  number  of  parameters  present  in  a  subset,  so  any  subset 
regression  model  contains  p-1  X  variables.  And  n  is  the  number  of 
observations . 

2  2 

The  Rp  Criteria.  Rp  is  based  on  the  coefficient  of  multiple 
2 

determination  R  ,  for  a  subset  of  size  p,  and  is  defined  as: 

Rp  =  1  -  (SSEp  /  SSTO)  (31) 

where 

SSEp  =  Error  sum  of  squares  for  a  parameter  subset  of  size  p. 

SSTO  =  Total  sum  of  squares  for  y. 

SSTO  is  equivalent  to  SSEj  which  is  the  regression  model  with  only  an 

intercept  term.  SSTO  remains  constant  for  each  subset  evaluated,  so  as 
2 

p  increases,  Rp  increases.  This  occurs  since  SSEp  can  not  increase  as 

2 

additional  variables  are  added  to  the  model.  Consequently,  Rp  reaches  a 

maximum  when  all  P-1  variables  enter  the  regression  model.  Therefore 

2 

the  intent  is  not  to  maximize  Rp,  but  to  find  a  point  where  adding 

2 

additional  variables  to  the  model  does  not  increase  Rp  significantly. 
'Clearly,  the  determination  of  where  diminishing  returns  sets  in  is  a 
judgmental  one'  (Neter,  et  al;  1983:422). 

2 

The  Ra  Criteria.  The  adjusted  coefficient  of  multiple 
2  2 

determination,  Ra,  is  very  similar  to  Rp,  and  is  defined  as: 


(32) 


Ra  =  1  -  ( (n-l/n-p) » (SSEp/SSTO) ) 

2 

However,  unlike  Rp,  this  criterion  'takes  the  number  of  parameters  in 

the  model  into  account  through  the  degrees  of  freedom*  (Neter,  et  al; 

2 

1983:424).  Therefore,  while  seeking  the  maximum  value  of  Ra,  it  is 
possible  for  this  value  to  decrease  as  p  increases  if  the  reduction  in 
SSGp  is  too  small  to  offset  the  loss  of  a  degree  of  freedom. 

The  Cp  Criteria.  The  Cp  criteria  is  based  on  minimizing  the 
total  mean  squared  error  of  the  n  fitted  values  for  each  of  the  various 
subset  regression  models.  Cp  is  defined  as: 

Cp  =  (SSEp  /  MSE  (X]  ,  .  .  .  ,  Xp_  2 ) )  -  (n  -  2*p)  (33) 

where 

MSE (Xi , . . ,Xp_ j)  a  mean  squared  error  of  the  model  with  all  P 
parameters . 

The  above  equation  assumes  that  the  model  which  includes  all  P 
parameters  provides  an  unbiased  estimate  of  the  variance.  In  the  event 
the  model  used  has  substantial  bias,  it  may  be  best  to  expand  the  set  of 
potential  variables  to  eliminate  the  bias. 

In  using  the  Cp  criterion,  identification  of  an  appropriate  subset 
of  X  variables  is  based  on:  (1)  A  small  value  of  Cp,  and  (2)  A  Cp  value 
near  p.  When  plotting  Cp  vs.  p,  models  with  little  bias  will  fall  near 
the  line  Cp  =  p,  and  models  with  significant  bias  will  be  substantially 
above  the  line. 

The  BCp  Criteria.  The  BCp ,  or  Best  Controls,  criteria  is  a 
new  criteria  developed  by  Bauer  and  Wilson  (1990).  Unless  otherwise 


18 


noted,  the  following  discussion  is  based  on  their  paper. 

Development  of  the  selection  criteria  is  presented  for  two  cases. 
The  first  case  is  where  the  covariance  matrix,  of  variables  and 
responses,  is  estimated.  The  second  case  is  where  the  covariance  matrix 
is  known. 

Before  proceeding,  the  ‘best*  subset  of  variables  is  defined  as 
that  subset  which  produces  a  confidence  region  of  minimum  expected 
volume.  Consequently,  the  selection  criteria  is  designed  to  identify 
which  subset  of  variables  will  produce  this  result. 

Nomenclature .  Let  Y  =  (Yj,...,Yp)’  denote  a  column 
vector  of  p  responses  generated  on  a  single  run  of  a  simulation  model 

ft 

whose  mean  response  Uy  =  E(Y)  is  to  be  estimated.  Furthermore,  let  C  = 
(Ci . Cq) ’  denote  a  column  vector  of  q  concomitant  control  variates 

ft 

with  known  mean  uc  =  E(C)  and  let  b  denote  a  fixed  (p  x  q)  matrix  of 
control  coefficients,  then  the  controlled  response 

Y (b)  =  Y  -  b(C  -  uc)  (34) 

is  an  unbiased  estimator  of  Uy  whose  dispersion  can  be  minimized  by  the 
appropriate  choice  of  b.  Let  Ey,  Ec,  and  Ey0  respectively  denote  the 
covariance  matrices  of  Y,  C,  and  between  Y  and  C;  then 


Ey  = 

Cov(Y)  =  E[(Y  -  Uy) (Y  -  Uy) ’ ), 

(35) 

Ec  * 

Cov(C)  >  E[(C  -  uc) (C  -  uc) 

(36) 

Eye 

=  Cov (Y ,C)  =  EUY  -  Uy)  (C  -  uc)  ’  ]. 

(37) 

And  then,  in  terms  of  these  quantities,  the  conditional  covariance 
matrix  of  Y  given  C  *  c  is 
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Ey  >  c  ~  Cov(Y!C  =  c)  =  Ey  ~  Eye  -  Eyc^c  ^yc  (38) 

q  q 

for  every  c  which  is  an  element  of  R’  (i.e.  cfcR  )  . 

Rubinstein  and  Marcus  (1985)  showed  that  the  generalized  variance  of 
Y(b)  is  minimized  by  the  optimal  matrix  of  control  coefficients 

M  •**  Al  1 

8  =  EycEc  (39) 

IV  M 

Typically,  EyC  is  unknown  so  8  must  be  estimated.  Let  k  denote  the 
number  of  independent  replications  of  the  simulation  to  be  performed; 

th 

and  for  j  =  1 . k,  let  (Yj.Cj)  denote  the  results  observed  on  the  j 

run.  Then,  in  terms  of  the  statistics 

k 


Y  =  ( 1/k)  E  Yj  , 
j  =  l 

k 

(40) 

Sy  =  ( 1/  (k- 1) )  E  (Yj  -  Y) (Yj  - 
j  =  l 

k 

Y)\ 

(41) 

C  =  (1/k)  E  Cj , 
j  =  l 

k 

(42) 

Sc  =  ( 1/ (k- 1 ) )  E  (Cj  -  C)  (Cj  - 
j=l 

and 

k 

C)\ 

(43) 

Syc  =  ( 1/ (k-1) )  E  (Yj  -  Y) (Cj 

J  =  1 

the  sample  analogue  of  8  is 

-  C)  ' ; 

(44) 

8  =  SycS^1 

(45) 

Thus  the  j*’*1  controlled  response  is  estimated  as 

Yj(8)  =  Yj 

-  8  (Cj  -  uc) 

for  J  =  1 . k;  and  the  overall  controlled  point  estimator  of  uy  is 
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Y(fi)  =  ( 1/k)  Z  Yj(fi)  =  Y  -  fi(C  -  uc)  (46) 

j  =  l 

Bauer  and  Wilson  (1990:4)  also  state  that  ‘In  large  scale 
simulation  experiments,  frequently  the  responses  and  the  controls  are 
jointly  normal  because  these  statistics  are  simultaneously  accumulated 
over  the  duration  of  each  run  and  thus  are  subject  to  a  central-limit 
effect;...  Thus  it  is  reasonable  to  assume  that  Y  and  C  jointly  possess 
a  multivariate  normal  distribution.*  Thus, 


Ey  Eyc 

Egy  Ec 


with  det(Ey)  f  0  and  det(Ec)  t  0,  where  uc  is  known  but  Ey,  EyC,  and 
possibly  Ec  are  unknown.  Applying  the  basic  results  of  Rao  (1967)  in  this 
situation,  Bauer  and  Wilson  then  compute  a  confidence  ellipsoid  for  uy 
as  follows.  Let 


T2(u)  =  [Y(fi) -u] ' [d’dEy |c3_ 1 [Y(fl) -u] 


for  every  u  which  is  an  element  of  R  ,  where 


(Ct  -  C) ’ 


(Ck  -  C)’ 


,  d  =  (1/k) lk  -  M(M’M) ” 1 (C  -  uc) , 


Ey  j  c  =  (k-1  /  k-q-l)(Sy  -  SycSc  S£c), 


and  let  lk  denote  a  k-dimensional  column  vector  of  ones.  Conditioned  on 

the  values  of  the  controls  (Cj:j  =  l,...,k)  observed  across  all  k  runs, 

2  "  2 

T  (Uy)  has  Hotelling’s  T  -distribution  with  k-q-1  degrees  of  freedom; 
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therfore  an  exact  unconditional  100(l-a)X  confidence  ellipsoid  for  uy  is 

M(q;k,p,a)  =  {u*RP : (T2 (u)  /  k-q-1)  *  (k-q-p  /  p) 

<  Fi-afp.k-q-p) ) ,  (51) 

where  Fi_a(p, k-q-p)  is  the  quantile  of  order  1-a  for  an  F-distribution 
with  p  and  k-q-p  degrees  of  freedom. 

Selection  Criteria  for  Estimated  Covariance  Matrix. 

Given  the  replication  count  k,  the  p-dimensional  estimand  uy,  and  the 
confidence  coefficient  ‘a‘,  the  goal  is  to  select  a  subset  of  controls 
from  a  set  of  q  control-variate  candidates  such  that  the  resulting 
controlled  confidence-region  estimator  for  uy  analogous  to  Eq  (51)  is, 
in  some  sense,  as  ‘small'  and  as  'stable'  as  possible.  Bauer  and  Wilson 
formulate  such  an  estimator  with  some  additional  notation.  Given  a 
nonnegative  integer  r  representing  the  number  of  control-variate 
candidates  currently  under  consideration,  let  u(q,r)  =  q ! / (r ! (q-r) ! ]  be 
the  number  of  distinct  control-variate  subsets  of  size  r.  Then,  for  r  = 
0 , . . . ,q  and  h  =  l,...,u(q,r)  let  I(h,r)  denote  the  hth  distinct  subset 

of  size  r  from  the  set  (1 . q) .  Furthermore,  on  the  run  of  the 

simulation  model,  let  Cj(h,r)  denote  the  r-dimensional  vector  of 
controls  corresponding  to  the  index-set  I(h,r) 

Cj(h.r)  =  [C(i i j) . C(irj) ] ’ ,  (52) 

where  ij  <  ...  <  ir  and  (ij,...,ir)  =  I(h,r). 

Similarly,  let  Ey!c(h,r) ,  fi(h,r),  Y[fi(h,r)],  d(h,r),  and  Ey;c(h,r) 
respectively  denote  the  analogues  of  Eqs  (38)  ,  (45) ,  (46)  ,  (49) ,  and 
(50)  when  the  control  vector  C(h,r)  defined  by  I(h,r)  is  used  to  compute 


the  controlled  estimator  of  Uy.  Corresponding  to  Eq  (48),  Bauer  and 
Wilson  then  derive 

T2  (u,h  ,r)  =  {YU(h,r)3  -  u) '  td  ’  (h  ,r)  d  (h  ,r)  Ey .  c  (h  ,r)  ]-1  * 

{Y[fi (h ,r)  ]  -  u);  (53) 

and  the  exact  lOO(l-a)*  confidence  ellipsoid  for  uy  analogous  to  Eq  (51) 
as 

M(h,r;k,p,a)  =  (u£RP:  (T2(u,h,r)  /  k-r-1)  *  (k-r-p  /  p) 

<  F1.a(p,k-r-p)).  (54) 

Then  the  size  of  the  confidence  region  is  calculated  by 

V(h,r;k,p,a)  =  ( j Ey . c(h ,r) | 1  /  (p/2) Q(p/2) )  * 

{ [d ’ (h,r)d(h,r) ] (pi«p* (k-r-1  /  k-r-p)  * 
F!..a(p,k-r-p)}p/2  (55) 

where  Q,  for  this  and  following  equations,  denotes  the  Gamma  function. 
Also  the  mean  volume  of  the  confidence  region,  Eq  (54),  is  given  by 

E(V(h,r;k,p,a)  ]  =  w(h,r;k,p,a)  ([G(k/2)  ]W2  /  G[  (k-p)/2J)  ,  (56) 

where 

w(h,r ;k,p,a)  =  ([ |Ey!e(h,r) |G(k/2) ]1/2  /  (p/2)G(p/2))  * 

C2*pi*p»F1.a(p, k-r-p)  /  k  (k-r-p) ]p/2 .  (57) 

And  the  mean  square  volume  of  the  confidence  region,  Eq  (54),  is 

E[V2(h,r;k,p,a)  ]  =  (w2(h,r;k,p,a)  /  G[(k-2p)/2])  « 

P 

fj  t (k-r-i) /  (k-r-2i)  ] .  (58) 

i  =  l 
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Although  alternative  expressions  are  available  for  the  mean  volume  and 
mean  square  volume,  Bauer  and  Wilson  believe  their  equations  (56)  and 
(58)  are  easier  to  use  from  both  a  mathematical  and  computational 
standpoint . 

In  the  case  of  a  univariate  response,  Nelson  (1989)  and  Schmeiser 
(1982)  define  the  standard  measures  of  confidence-interval  stability  as 
the  standard  deviation  (SD)  and  coefficient  of  variation  (CV)  of  the 
confidence-interval  half-length;  and  in  the  case  of  a  multivariate 
response,  the  corresponding  stability  measures ,  based  on  the  confidence- 
region  volume  of  Eq  (55)  ,  are 

SD[ V(h ,r ; k , p , a) ]  =  w(h , r ;k ,p ,a)  {(1  /  G[(k-2p)/2])  « 

fj  [(k-r-i)/(k-r-2iJ j  -  (G(k/2)  /  G2 C (k-p) /2 )) ) 1/2  (59) 

i  =  1 

and 

CV[ V (h ,r ; k ,p ,a) ]  =  ( (G2[ (k-p) /2]  /  G(k/2) G[ (k-2p) /2 ] )  * 

/7  t(k-r-i)/(k-r-2i)  ]  -  1}1/2.  (60) 

i  =  l 

Similar  to  Schmeiser’s  (1982)  conclusions  about  the  performance  of 
univariate  confidence  intervals,  Bauer  and  Wilson  (1990)  observed  that  a 
confidence-region  estimator  M(h,r;k,p,a)  with  large  values  of  (59)  or 
(60)  will  give  false  signals  about  the  intrinsic  precision  of  the 

A 

Ar 

associated  point  estimator  Y[6(h,r)3  in  a  large  percentage  of 

applications.  Bauer  and  Wilson  then  went  on  to  state 

Thus  it  seems  reasonable  to  select  a  control  vector  C(h,r)  that 
yields  a  small  value  for  (59)  or  (60).  However,  it  would  be 
undesirable  to  base  a  control-variate  selection  criterion 
exclusively  on  the  principle  of  minimizing  (59)  or  (60)  --  this 
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principle  fails  to  exclude  confidence-region  estimators  that 
achieve  smaller  values  of  the  standard  deviation  or  coefficient  of 
variation  of  the  volume  simply  by  increasing  the  mean  volume.  On 
the  other  hand,  it  is  also  undesirable  to  base  a  control-variate 
selection  criterion  exclusively  on  the  principle  of  minimizing  the 
mean  volume  without  attempting  simultaneously  to  reduce  or  at  least 
bound  the  standard  deviation  or  coefficient  of  variation  of  the 
volume.  (1990:6) 

Bauer  and  Wilson  wanted  to  ensure  that  the  delivered  confidence- 

region  estimator  Eq  (54)  is  both  small  and  stable,  so  they  proposed  a 

control-variate  selection  criterion  based  on  the  principle  of  minimizing 

2 

the  mean  square  volume  given  by  Eq  (58).  Since  E[V  (h,r;k,p,a)]  > 

E^[V(h,r;k,p,a) ]  >  0,  it  is  clear  that  their  selection  criterion  will 

tend  to  reduce  the  mean  volume  at  least  indirectly;  and  in  comparison  to 

a  selection  criterion  based  on  minimization  of  the  mean  volume,  Eq  (56), 

the  selection  criterion  will  yield  smaller  values  of  the  standard 

deviation  Eq  (59)  and  the  coefficient  of  variation  Eq  (60)  of  the 

confidence-interval  volume;  unless  both  procedures  select  exactly  the 

same  control  variates  with  probability  one.  In  the  event  both 

procedures  do  select  the  same  control  variates,  the  two  selection 

criteria  yield  identical  results.  Therefore  their  strategy  of  basing 

the  criterion  on  minimizing  the  mean  square  volume  of  the  delivered 

confidence-region  estimator  offers  many  of  the  advantages  of  selection 

criteria  based  on  minimizing  the  mean,  standard  deviation,  or 

coefficient  of  variation  of  the  delivered  volume  without  some  of  the 

potential  disadvantages  of  these  latter  selection  criteria. 

To  implement  the  proposed  selection  criterion  in  practice,  it  is 

2 

necessary  to  minimize  the  mean  square  volume  E[V  (h,r;k,p,a)]  as  a 

function  of  h  and  r,  where  r  =  0,...,q  and  h  =  1 . u(q,r).  Since 

IEy;c(h,r);  is  generally  unknown,  this  quantity  is  replaced  by  the 
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unbiased  estimator 


~  P 

|Ey!c(h,r) |  /7  I (k-r- 1) / (k-r-i) ]  (61) 

i  =  l 

to  obtain  the  expression  that  must  be  minimized  in  selecting  the  final 
subset  of  control  variates 


MIN  C ( I Ey ; c ( h , r > I G(k/2)  /  [ (p/2) G(p/2) ]2G[ (k-2p) /21) 

»  [ (2«pi*p»Fi_a(p,k-r-p)  /  k(k-r-p)]^ 

P 

*  77  [ (k-r-l) / (k-r-i)  ])  (62) 

i  =  1 

ft  * 

subject  to  the  constraints  of  0<r<q  and  l<h<u(q,r).  Let  r  and  h 
denote  the  optimal  values  of  r  and  h  in  Eq  (62) .  The  delivered  point 

~  ~  *  ft 

and  confidence-region  estimators  of  uy  are  given  by  Y[fi(h  ,r  )]  and 

*  ft 

M(h  ,r  ;k,p,a),  respectively.  Thus  Eq  (62)  gives  the  selection  criteria 
for  the  case  where  the  covariance  matrix  is  unknown. 

Selection  Criteria  for  Known  Covariance  Matrix.  Often, 
situations  arise  in  discrete  event  simulation  where  the  covariance 
matrix  of  some  set  of  control  variates  is  known  analytically  or  can  be 
readily  evaluated  by  numerical  methods.  In  this  situation  an 
alternative  to  the  estimator  &,  for  the  unknown  covariance  matrix  case, 
for  the  optimal  control  coefficient  vector  fi  can  be  obtained  by 
replacing  Sc,  in  Eq  (45),  with  Ec  to  obtain 


sycEc 


(63) 


In  this  case  the  controlled  point  estimator  of  uy  has  the  form 

t  *. 

Y(B)  =  Y  -  fi(C  -  uc). 


(64) 
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t 

Under  the  assumption  of  joint  normality,  Bauer  (1987)  proved  that  Y(fi) 
is  an  unbiased  estimator  of  uy  with  covariance  matrix 

•  • 

E  =  Cov[  Y (fl) ]  =  (k-2  /  k (k- 1 ) ) Ey ; c  ♦  (q+1  /  k(k-l))Ey.  (65) 

To  derive  an  approximate  100{l-a)X  confidence  region  for  uy 
* 

centered  at  Y(fi),  Bauer  and  Wilson  began  with  an  unbiased  estimator  of 
obtained  by  replacing  the  unknown  covariances  on  the  right-hand  side  of 
Eq  (65)  with  the  corresponding  sample  covariances 

S  =  (k-2  /  k(k-l) )Ey;c  +  (q+1  /  k (k- 1) ) Ey .  (66) 

Provided  that  k  >>  q,  Eq  (66)  implies  that  S  is  approximately 

t 

independent  of  Y(fl)  and  possesses  the  p-dimensional  central  Wishart 
distribution  with  k-q-1  degrees  of  freedom  on  the  covariance  matrix 
Ey • c  .  Then 

T2 (u)  =  (Y(fi)  -  u)’S_1[Y(fi)  -  u]  (67) 

~  p  *2  ~ 

for  every  u  which  is  an  element  of  R  .  Thus  for  large  k,  T  (uy)  has 

2 

an  approximate  Hotelling’s  T  -distribution  with  k-q-1  degrees  of  freedom; 
and  in  this  case  an  approximate  100(l-a)%  confidence  ellipsoid  for  uy  is 
given  by 

M(q;k,p,a)  =  (u4Rp: (T2 (u)  /  k-q-1) (k-q-p  /  p) 

<  Fj.a(p, k-q-p)).  (68) 

Progress  to  this  point  parallels  the  development  of  equations  (54) 
through  (62)  that  apply  to  the  situation  where  the  covariance  matrix  is 
unknown.  Bauer  and  Wilson  then  sought  a  control  vector  C(h,r)  which 
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minimized  the  mean  sqaure  volume  of  the  confidence  ellipsoid.  For  r  = 

0 . q  and  h  =  1 . u(q,r),  let  fl(h,r),  Ytfl(h.r)],  E(h,r),  and  S(h,r) 

respectively  denote  the  analogues  of  Eqs  i63)  ,  (64) ,  (65) ,  and  (66)  when 
the  control  vector  C(h,r)  defined  by  the  index-set  I(h,r)  is  used  to 
compute  the  controlled  estimator  of  Uy.  Equation  (67),  then  becomes 

T2(u,h,r)  =  {Y[fi(h,r)3  -  u) '  [S  (h  ,r)  f 1  {YU (h  ,r) )  -  u)  (69) 

for  every  u  which  is  an  element  of  RP;  and  the  approximate  100(l-a)X 
confidence  ellipsoid  for  Uy  analogous  to  Eq  (68)  is 

M(h,r;k,p,a)  =  {u£RP: (T2 (u,h ,r)  /  k-r-l)(k-r-p  /  p) 

<  F)-a(P.k-r-p)) .  (70) 

Now  the  confidence  region  of  Eq  (70)  has  volume 

V(h ,r ;k ,p ,a)  =  (|S(h,r)(1/2  /  Cp/2) G (p/2) )  » 

[pi*p* (k-r- 1)  /  k-r-pJFj-alp.k-r-p) ]p/2  (71) 

• 

Then,  assuming  S(h,r)  approximately  possesses  the  p-dimensional  Wishart 
distribution  with  k-r-1  degrees  of  freedom  on  the  covariance  matrix 
E(h,r),  Bauer  and  Wilson  derived  the  following  approximate  expression 
for  the  mean  square  volume 

E[V2 (h ,r ;k ,p,a)  2  (|E(h,r)|  /  ( (p/2) G(p/2) ]2)  * 

p 

Cpi»p  /  k-r-p)F2_a(p,k-r-p) ]p/7  (k-r-i).  (72) 

i=  1 

To  implement  their  proposed  selection  criterion  for  this  new  linear 
control-variates  estimation  procedure,  it  is  necessary  to  minimize  Eq 
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(72)  as  a  function  of  h  and  r  in  a  manner  similar  to  that  for  Eq  (62)  . 
Since  !S(h,r)!  is  generally  unknown,  this  quantity  is  replaced  by  the 
unbiased  estimator 

.  P 

|S(h,r) j  /7  C (k-r- 1) / (k-r-i )  ]  (73) 

i  =  1 

as  in  Eq  (61)  to  obtain  the  expression  that  must  be  minimized  in 
selecting  the  final  subset  of  control  variates 

MIN  ( ( | S (h ,r) I  /  E (p/2)Q(p/2) )2)  » 

tpi»p«(k-r-l)*Fi-a(p,k-r-p)  /  k-r-p]p)  (74) 

#  *  *  ft 

subject  to  the  constraints  of  0<r<q  and  l<h<u(q,r).  If  r  and  h 
denote  the  optimal  values  of  r  and  h  in  Eq  (74)  then  the  delivered  point 

~  -  *  *t  #» 

and  confidence-region  estimators  of  uy  are  YEfi(h  ,r  ))  and 

*  *  *  *  * 

M(h  ,r  ;k,p,a),  respectively.  Thus  Eq  (74)  gives  the  selection 
criteria  for  the  case  where  the  covariance  matrix  is  unknown. 

Selection  Procedures.  There  are  a  variety  of  procedures  available 
for  selection  of  the  best  or  near-best  subset (s)  of  variables.  The 
procedures  presented  here  are:  enumerated  subsets,  stepwise  regression, 
forward  selection  and  backward  elimination,  and  regression  by  leaps  and 
bounds . 

Each  procedure  has  its  own  advantages  and  disadvantages.  However, 

the  primary  rationale  for  using  selection  procedures  other  than  an 

enumerated  subsets,  or  all-possible  subsets,  approach  is  to  reduce  the 

amount  of  computations  required.  Even  for  a  moderate  number  of 

p 

variables  the  number  of  subsets  to  evaluate  becomes  2  -  1,  where  P  is 

the  number  of  variables  being  considered  for  the  model.  In  addition, 
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only  the  enumerated  subsets  approach  will  assure  that  the  'best' 
subset (s)  of  variables  will  be  selected,  based  on  the  selection  criteria 
used . 

Enumerated  Subsets.  The  enumerated  subsets,  or  all-possible 
subsets,  approach  is  based  on  the  following  algorithm. 

1)  A  regression  model  with  no  X  variables  (i.e.  +  e^)  is 

evaluated  using  the  selected  criteria. 

2)  A  series  of  regression  models  including  each  variable,  individually, 
are  evaluated  using  the  selected  criteria. 

3)  A  series  of  regression  models  including  each  possible  pair  of 
variables  are  evaluated.  This  process  continues,  increasing  the 
number  of  variables  in  the  model  one  at  a  time,  until  a  model 
incorporating  all  variables  is  reached.  Again,  each  model  is 
evaluated  using  the  selected  criteria. 

2  2 

4)  Based  on  the  criteria  employed  (i.e.  Rp,  Ra,  Cp,  BCp,  etc.),  the 
best,  or  k  best,  subset (s)  are  selected. 

Stepwise  Regression.  The  stepwise  regression  procedure  is  the 

most  common  search  method  used  when  the  number  of  variables  involved 

make  an  enumerated  subsets  approach  infeasible. 

Essentially,  this  search  method  develops  a  sequence  of  regression 
models,  at  each  step  adding  or  deleting  an  X  variable.  The 
criterion  for  adding  or  deleting  an  X  variable  can  be  stated 
equivalently  in  terms  of  error  sum  of  squares  reduction, 
coefficient  of  partial  correlation,  or  F  statistic.  (Neter,  et  al; 
1983:430) 

ft 

For  the  search  algorithm  which  follows,  the  F  statistic  is  used 
for  illustration  of  the  concepts  involved. 

1)  A  regression  model  is  fitted  for  each  of  the  X  variables.  For  each 

* 

regression  model  the  F  statistic  is  obtained  as  follows: 
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(75) 


fJ  =  MSR(Xk)  /  MSE(Xk) 

ft 

The  X  variable  with  the  largest  F  value  is  selected  to  enter  the 

* 

model.  Provided  the  F  value  exceeds  a  predetermined  level  required 
to  enter  the  model.  If  all  of  the  F*  values  are  below  the  threshold 
level,  the  search  ends. 

2)  If  an  X  variable,  say  X^ ,  enters  the  model,  then  all  regression 
models  with  two  variables  are  fitted;  where  Xj  is  one  of  each  pair. 
For  each  regression  model,  the  partial  F  test  statistic  (Neter,  et 
al;  1983:289)  is  obtained. 

Fj  =  MSR(Xk:Xi)  /  MSEfXi.Xjt)  =  (bk  /  s(bk))2  (76) 

Where  bk  and  s(bk)  are  the  estimated  regression  coefficient  of 
variable  k  and  the  estimated  variance  of  the  estimated  regression 
coefficient  of  variable  k,  respectively.  Again,  the  X  variable  with 

ft 

the  largest  F  value  enters  the  model,  if  it  exceeds  the  threshold 
level.  Otherwise  the  procedure  ends. 

3)  When  more  then  one  variable  enters  the  model,  it  is  then  determined 

ft 

if  any  of  the  variables  in  the  model  should  be  dropped.  The  F 
values  are  derived  as  follows: 

Fk  =  MSR(Xk!Xi . Xj )  /  MSE(Xk,Xi . Xj )  (77) 

Where  Xk  is  the  variable  being  tested  for  possible  elimination  from 
the  model  and  X^,...,Xj  are  the  other  variables  in  the  model.  The  X 
variable  with  the  smallest  F*  value  is  selected  to  exit  the  model. 

ft 

Provided  the  F  value  falls  below  a  predetermined  level  required  to 
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exit  the  model.  If  all  of  the  F*  values  are  above  the  threshold 
level,  all  variables  remain  in  the  model. 

4)  Steps  (2)  and  (3)  are  repeated  until  no  further  X  variables  can  meet 
the  threshold  levels  to  enter  or  exit  the  model.  It  should  be  noted 
that  as  the  number  of  variables  in  the  model  increases,  the  size  of 
the  subsets  evaluated  in  step  (2)  increases.  Where  each  subset 
evaluated  includes  the  variables  currently  in  the  model,  plus  one  of 
the  variables  not  in  the  model. 

Forward  Selection  and  Backward  Elimination.  Both  forward 
selection  and  backward  elimination  procedures  are  simplified  variations 
of  the  stepwise  regression  procedure.  The  forward  selection  procedure 
differs  from  the  stepwise  regression  procedure  by  not  testing  a 
variable,  once  it  has  entered  the  model,  for  possible  elimination  from 
the  model . 

The  backward  elimination  procedure  is  the  opposite  of  the  forward 

selection  procedure.  This  procedure  begins  with  the  model  containing 

* 

all  the  X  variables.  Then  the  F  value  for  each  variable  is  calculated 
and  the  variable  with  the  smallest  value  identified.  If  this  value  is 
less  than  the  threshold  level,  it  is  dropped  from  the  model.  This 
process  continues  until  no  further  variables  can  be  dropped  from  the 
model . 

From  this  overview  of  control  variates  theory  and  methods  for 
selecting  significant  variables,  a  methodology  was  developed  for 
proceeding  with  the  work  involved.  This  methodology  is  presented  in  the 
following  chapter. 
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III.  METHODOLOGY 


The  overall  objective  of  this  thesis  was  to  develop  software  to 
assist  in  identifying  the  significant  control  variables  in  a  simulation 
model.  The  software  was  to  incorporate  a  newly  developed  selection 
criteria  in  conjunction  with  several  common  selection  procedures.  After 
this  was  accomplished,  the  objective  moved  to  applying  the  software  to  a 
simulation  model  and  evaluating  the  results. 

Selection  and  Preparation  of  Computer  Code 

The  first  step  in  proceeding  was  to  select  and  prepare  the 
necessary  computer  code/software.  The  software  specifically  is  the 
Variable  Subset  Selection  Program,  the  simulation  model  which  provided 
data  to  test  it  on,  and  software  for  processing  the  data  collected  from 
the  simulation  model. 

Variable  Subset  Selection  Program.  The  Variable  Subset  Selection 
Program  was  developed  from  a  previously  written  program.  Before  the 
software  was  ready  for  use,  extensive  revision  and  expansion  of  the  code 
was  performed.  The  primary  goals  achieved  in  revising  the  software 
dealt  with  increasing  the  ease  of  use  of  the  software,  and  permitting 
either  manual  or  external  file  data  entry.  Prior  to  these  revisions, 
the  program  data  had  to  be  coded  directly  into  the  software  before  it 
could  be  compiled  and  executed. 

Expansion  of  the  software  dealt  with  the  addition  of  a  capability 
to  perform  a  stepwise  regression  procedure  in  conjunction  with  the  new 
selection  criteria.  This  was  in  addition  to  the  enumerated  subsets 
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procedure  originally  incorporated. 

Data  Generation  Model.  Additionally,  a  simulation  model  was 
selected  to  generate  data  for  testing  the  variable  subset  selection 
software.  It  was  decided  to  select  a  model  which  simulates  the 
operation  of  a  Local  Area  Network  (LAN)  and  the  interaction  of  the 
system  peripherals.  The  network  model,  on  which  this  simulation  is 
based,  is  commonly  found  in  simulation  journals  and  serves  as  a 
practical  benchmark  for  testing  purposes. 

The  LAN  simulation  model  consists  of  seven  nodes  and  is  illustrated 
in  Figure  3.1.  Node  1  represents  the  terminals  connected  to  the  LAN 
system.  Commands  to  the  LAN  system  are  generated  according  to  an 
exponential  distribution.  The  commands  are  received  at  node  2,  which 
serves  as  a  delay  buffer  and  simulates  the  possibility  of  all  system 
peripherals  being  busy  and  unable  to  accept  new  commands.  Next  the 
commands  move  to  node  3,  which  routes  the  command  to  a  system 
peripheral.  The  routing  of  a  command  from  a  node  to  any  other  node  in 
the  system  is  based  on  probabilities  contained  in  a  probability  matrix. 

A  command  routed  from  node  3  can  go  to  node  1,  4,  5,  6,  or  7.  The  queue 

capacity  of  nodes  4,  5,  6,  and  7  is  infinite  and  the  queue  discipline  is 

First-In-First-Out  (FIFO).  At  nodes  4,  5,  6,  and  7  the  service  times 
are  distributed  exponentially.  Once  a  command  is  through  being  serviced 
by  a  peripheral  unit,  it  returns  to  node  3,  where  it  is  routed  to 
another  peripheral  for  further  processing  or  sent  to  node  1.  When  a 
command  returns  to  node  1,  the  time  it  took  to  get  through  the  system  is 

noted  and  the  command  terminated.  The  SLAM  and  associated  FORTRAN  code 

for  this  simulation  model,  along  with  the  following  post-processing 
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code,  are  provided  in  Appendix  C. 


The  data  collected  using  this  model  was  related  to  the  number  of 
commands  through  the  system,  average  utilization  and  total  service  time 
at  each  system  node,  number  of  departures  from  a  node,  and  number  of 
departures  from  node  3  to  the  other  system  nodes.  This  data  was  then 
processed  for  use  by  the  Variable  Subset  Selection  Program. 

Data  Post-Processing  Software.  After  collection  of  output  data 
from  the  Data  Generation  Model,  data  post-processing  was  performed.  The 
purpose  of  this  post-processing  was  to  convert  the  control  variate  data 
into  the  form  of  work  and  routing  variables.  These  types  of  variables 
are  desirable  because  of  the  characteristics  they  possess.  The  desired 
characteristics  are  a  sample  mean  of  zero  and  sample  variance  of  one. 

The  post-processing  software  performed  this  conversion  according  to 
the  following  equations 

_  1  /O  f  ^  ,  t) 

Xk(t)  =  (f(k,t)  1  /  WfcffU)))  Z  (Uj  Ck)  -  uk)  /  sk  .  (78) 

j  =  l 


where 


f(k,t)  -  number  of  service  times  that  are  finished  at  node  k 
during  time  period  (0,t) 

wk  s  relative  frequency  with  which  a  customer  visits  node  k 
Uk(j)  -  the  j-th  service  time  at  node  k 
uk  *  E[Uk(j) ] 
sk  *  Var[Uk ( j ) ] 


and 


H(t) 

xk  *  L  t(uk(j)  -  pk (*) ) 

j  =  i 


(N(t) ( 1  -  pk(»))pk(»))"1/2]  (79) 


for  J  8  1,...,S  and  N(t)  >  0 
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where 


Pfcf*)  =  probability  of  transition  from  node  3  to  node  k 
UjfU)  =  flag  whether  or  not  transition  j  was  to  node  k.  If  so, 
then  ( j ) = 1 ,  otherwise  it  equals  zero. 

N(t)  =  total  number  of  transitions  from  node  3  up  to  time  t 
S  =  total  number  of  nodes  in  the  LAN  model. 

Equation  (78)  applied  to  the  work  variables  and  Eq  (79)  to  the  routing 
variables.  Additional  information  on  work  variables  can  be  found  in 
Lavenberg,  et  al  (1982),  and  more  information  concerning  routing 
variables  in  Bauer  (1987). 

The  Response  Variables 

Next,  a  response  vector  was  chosen  for  the  analysis.  The  response 
chosen  to  form  the  response  vector  was  the  time  it  takes  a  command  to 
pass  through  the  LAN  system.  The  time  through  the  system  begins  when 
the  command  is  issued  at  a  terminal  and  ends  when  a  command  returns  to 
node  1,  informing  it  that  the  command  has  been  executed. 

The  Control  Variables 

Next,  the  control  variables  were  chosen,  to  provide  a  pool  of 
variables  for  the  Variable  Subset  Selection  Program  to  select  from. 

The  control  variables  chosen  were  the  total  service  times  accumulated  at 
nodes  1,  3,  4,  5,  6,  and  7;  and  the  number  of  departures  from  the  CPU 
(node  3)  to  nodes  1,  4,  5,  6,  and  7  respectively.  This  provided  a  pool 
of  eleven  variables. 

The  variables  analogous  to  total  service  time  at  node  2,  and  number 
of  departures  from  the  CPU  to  node  2  and  node  3  were  not  used.  These 
variables  were  infeasible  since  their  corresponding  values  were  always 
zero.  This  resulted  from  the  probability  matrix  associated  with  moving 
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from  one  node  to  another. 

Tit*  selection  of  these  variables  also  provided  a  mix  of  work  and 
routing  variables  for  the  Variable  Subset  Selection  Program  to  select 
from.  The  service  time  variables  became  work  variables,  and  the 
departure  variables  became  routing  variables. 

The  Experimental  Procedure 

There  were  two  distinct  phases  to  the  experimental  procedure 
employed  in  this  thesis  effort.  The  first  phase  involved  testing  the 
Variable  Subset  Selection  Program  using  data  generated  with  known 
covariances  between  the  control  variables  and  known  significant  control 
variables.  The  second  phase  dealt  with  testing  the  Variable  Subset 
Selection  Program  on  output  from  an  untested  simulation  model  and 
evaluating  this  data/output  using  the  VSSP  and  SAS  software  on  the  AFIT 
VAX  system. 

Evaluation  Using  Known  Data.  The  known  data  was  previously  derived 
by  Bauer  and  Wilson  (1990)  in  testing  of  the  original  version  of  the 
selection  software.  A  series  of  Monte  Carlo  experiments  were  performed 
to  derive  the  data.  Bauer  and  Wilson  choose  to  use  five  control  variate 
candidates  and  two  responses.  Two  of  the  five  control  variates  were 
constructed  to  be  uncorrelated  to  the  responses;  therefore  these  control 
variates  acted  as  decoys.  Bauer  and  Wilson  extended  Eq  (47)  by 
partitioning  the  vector  of  control  variates  as  C  =  [X’  W’ ] ,  where  X’  = 

[ C i ,  C2,  C3]  consists  of  the  three  effective  control  variates  and  W’  * 
[C4,  Cg]  consists  of  the  decoy  control  variates.  This  resulted  in  the 
overall  stochastic  structure  of 
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(80) 


where  0  is  an  appropriately  dimensioned  matrix  of  zeros.  Three  versions 
of  the  covariance  structure  in  Eq  (78)  were  constructed.  In  all  three 
cases 


1.0  0.3 

1  0  0 

1  0 

= 

.  Ex  * 

0  1  0 

,  and  Ew  = 

0.3  1.0 

1 

0  0  1 

0  1 

(81) 


The  cross-correlation  matrix  between  the  responses  and  effective  control 
variates  was  of  the  form 


a  0  c 
0  b  c 


(82) 


The  choices  of  a,  b,  and  c  for  each  case  are  summarized  in  Table  3.1. 
Table  3.1  Cross-correlations  in  Eq  (82) 


Parameter 

I 

Case 

II 

III 

a 

0.8 

0.5 

0.2 

b 

0.5 

0.5 

0.2 

c 

0.5 

0.5 

0.2 

Bauer  and  Wilson  (1990) 


One  thousand  sets  of  data  were  then  derived  from  the  normal 
distribution,  Eq(80)  for  each  case.  Bauer  and  Wilson  then  evaluated  the 
resultant  data  sets  at  k  s  number  of  replications  *  10  and  m  =  number 
of  meta-experiments  =  100,  and  k  =  20  and  m  =  50  to  obtain  coverage  and 
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volume  reduction  percentages  for  the  delivered  confidence  regions. 

Once  the  Variable  Subset  Selection  Program  was  ready  and  the 
data  sets  obtained,  evaluation  was  performed.  Each  data  set  was 
evaluated  for  the  same  k  and  m  values,  used  by  Bauer  and  Wilson,  using 
both  the  enumerated  subsets  and  stepwise  procedures.  In  addition,  for 
this  thesis,  further  evaluation  of  the  original  data  was  performed  at  k 
=  50  and  m  =  20,  and  k  =  100G  and  m  =  1.  Again,  evaluation  was 
performed  using  both  enumerated  subsets  and  stepwise  procedures. 

Evaluation  Using  Untested  Simulation  Model.  This  section  outlines 
the  steps  followed  in  generating  the  data,  using  the  untested  simulation 
model,  for  this  thesis  project.  In  generating  the  data/output  from  the 
untested  simulation  model,  1000  runs  with  different  seeds  for  the  random 
number  generators  were  performed.  The  data  generated  from  each  run  of 
the  simulation  was  placed  into  an  output  file  called  'DQM.OP'.  The 
simulation  model  was  run  on  the  AFIT  VAX  system. 

After  the  data/output  was  collected  from  the  simulation  model,  the 
data  was  converted  into  the  appropriate  work  and  routing  variables  using 
the  data  post-processing  software.  As  the  work  and  routing  variable 
data  was  created,  it  was  placed  in  an  output  file  called  'VSSP1.IN*. 
Also,  a  datafile  called  ‘VSSP0.IN’  was  prepared  which  contained  any 
other  information  required  by  the  Variable  Subset  Selection  Program. 
'VSSPO.ir  and  ‘VSSP1.IN*  then  became  data  input  files  for  the  Variable 
Subset  Selection  Program. 

Next,  the  Variable  Subset  Selection  Program  was  run  using  data 
contained  in  the  data  input  files  ’VSSP0.IN*  and  ’VSSP1.IN'.  The  output 
of  the  software  was  then  placed  in  a  file  called  ‘VSSP.OUT*  and 
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contained  the  selected  subset (s)  of  the  control  variables,  deemed  to  be 
significant.  After  this  point,  the  data  evaluation  phase  occurred. 

Data  Evaluation  Procedure 

After  the  data  had  been  generated  and  collected  it  then  had  to  be 
evaluated.  For  the  results  obtained  using  the  known  data,  the 
evaluation  primarily  consisted  of  comparing  the  results  originally 
derived  by  Bauer  and  Wilson  were  compared  to  those  derived  using  the 
VSSP.  In  addition,  for  results  obtained  from  data  runs  beyond  those 
performed  by  Bauer  and  Wilson,  the  results  were  examined  to  see  what 
would  happen  to  the  results  and  if  any  noticeable  trends  would  develop. 

In  regards  to  results  obtained  using  the  untested  simulation  model. 
The  primary  evaluation  was  based  on  comparison  of  the  subsets  selected 
as  a  result  of  each  selection  procedure  employed,  the  BCp  criterion 
value  for  each  resulting  variable  subset,  and  the  estimated  coverage  and 
variance  reduction  associated  with  the  control  variables  selected. 

In  addition,  further  comparison  was  performed  against  results 
derived  using  the  SAS  enumerated  subsets  and  stepwise  regression 
procedures.  This  primarily  served  the  purpose  of  validating  the  results 
derived  with  the  BCp  criterion. 

The  results  derived  from  all  these  procedures  are  summarized  and 
discussed  in  the  following  chapter. 


41 


IV.  RESULTS  AMD  DISCUSSION 


The  results  obtained  from  this  effort  are  presented  and  discussed 
according  to  phase  of  the  experimental  procedure  in  which  they  were 
derived.  The  phases  dealt  with  evaluation  of  the  data/output  with  known 
characteristics,  and  evaluation  of  data/output  from  the  untested 
simulation  model. 

Results  From  Known  Data 

The  results  of  the  Variable  Subset  Selection  Program  runs,  using 
the  data/output  with  known  characteristics,  are  presented  in  Tables  4.1, 
4.2,  4.3,  and  4.4  for  the  various  combinations  of  k  and  m  tested. 

Examination  of  the  results,  using  the  data  originally  derived  by 
Bauer  and  Wilson  (1990),  reveals  several  items  of  interest.  First,  as 
the  number  of  replications  increased,  this  had  significant  effects  on 
the  results.  The  percentage  coverage  and  percentage  confidence  volume 
reduction  figures  stabilized  as  the  number  of  replications  increased. 
Also,  for  each  data  case,  these  figures  become  uniform,  regardless  of 
the  evaluation  method  used  and  whether  the  covariances  between  the 
control  variables  was  estimated  or  known.  Apparently,  some  asymptotic 
point  was  reached  where  increasing  the  number  of  replications  per  meta- 
experiment  ceased  to  make  a  difference. 

It  was  also  noted  that  for  the  initial  results,  for  ks10  and  m=100, 
the  evaluations  using  the  known  covariance  matrix  had  better  coverage. 
However,  this  advantage  quickly  disappeared  as  the  number  of 
replications  increased. 
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Table  4.1  VSSP  results  for  k  =  10,  m  =  100,  and  a  =  0.10 


Evaluation 

Method 

Selection 

Criterion 

X 

I 

Cover* 
ror  Caj 
II 

ige 

je 

III 

X  Vo] 
I 
I 

Bedu< 
ror  Cat 
II 

:tion 

?e 

III 

Enumerated 

Subsets 

Eq  (62) 

78 

76 

69 

96 

83 

74 

Eq  (74) 

94 

88 

81 

83 

75 

72 

Stepwise 

Eq  (62) 

77 

79 

71 

96 

83 

74 

Eq  (74) 

93 

88 

81 

83 

75 

71 

Table  4.2  VSSP  results  for  k  =  20,  m  =  50,  and  a  =  0.10 


Evaluation 

Method 

Selection 

Criterion 

X 

l 

Coveri 
ror  Cai 
II 

ige 

;e 

III 

X  VoJ 
I 
I 

Reduc 
'or  Cai 
II 

:tion 

se 

III 

Enumerated 

Subsets 

Stepwise 

Eq  (62) 

84 

80 

82 

95 

76 

65 

Eq  (74) 

Eq  (62) 

82 

84 

84 

80 

II 

II 

QO  II  00 

M  II  dfc 

II 

88 

95 

74 

76 

64 

65 

Eq  (74) 

82 

84 

84 

88 

74 

64 

Next  it  was  noticed  that  for  the  first  two  sets  of  results  (k=10  & 
m=100,  and  k=20  &  m=50) ,  the  results  between  the  known  and  estimated 
covariance  runs,  appears  to  be  converging.  When  the  estimated 
covariances  were  used,  coverage  started  out  low  and  vice  versa  when  the 
known  covariances  were  used.  However,  later  sets  of  results  rebuked 
this  observation.  It  is  not  yet  understood  why  this  occurred. 

Finally,  it  was  found  that  the  volume  reduction  achieved  was 
greater  for  data  cases  with  greater  covariances  between  the  control 
variables  and  responses.  This  becomes  more  pronounced  as  the  number  of 
replications  increase,  but  is  still  readily  apparent  even  for  a  low 
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number  of  replications.  Conversely,  the  covariance  structure  of  the 


data/output  had  little,  if  any,  impact  on  the  coverage  figures.  This 
result  makes  intuitive  sense  since  as  covariances  increase,  the 
data/output  becomes  more  tightly  grouped;  thus  more  data  would  fall 
within  the  calculated  confidence  volume. 


Table  4.3  VSSP  results  for  k  =  50,  m  *  20,  and  a  =  0.10 


Evaluation 

Method 

Selection 

Criterion 

X 

I 

I 

Cover* 
'or  Ca£ 
II 

ige 

le 

III 

X  Vo) 

I 

Reduc 
'or  Cas 
II 

:tion 

se 

III 

Enumerated 

Subsets 

Eq  (62) 

75 

75 

75 

95 

76 

61 

Eq  (74) 

75 

70 

75 

92 

74 

60 

Stepwise 

Eq  (62) 

75 

75 

75 

95 

76 

61 

Eq  (74) 

75 

II 

II 

II  -J 

II  o 

II 

75 

92 

74 

60 

Table  4.4  VSSP  results  for  k  =  1000,  m  =  1,  and  a  =  0.10 


Evaluation 

Method 

Selection 

Criterion 

X 

l 

Cover* 
ror  Cat 
II 

ige 

se 

III 

X  Vo] 

I 

Reduc 
ror  Cas 
II 

:tion 

se 

III 

Enumerated 

Subsets 

Eq  (62) 

100 

100 

100 

95 

76 

59 

Eq  (74) 

100 

100 

100 

95 

76 

59 

Stepwise 

Eq  (62) 

100 

100 

100 

95 

76 

59 

Eq  (74) 

100 

100 

100 

95 

76 

59 

Additionally,  the  variable  subsets  selected  by  each  evaluation 
method  were  compared  for  the  various  combinations  of  k  and  m  tested. 
The  purpose  of  this  was  to  identify  if  any  trends  developed.  The 
results  are  presented  in  Table  4.5. 
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Table  4.5  Percentage  Number  o f  Different  Variable  Subsets 

Selected;  Comparison  Between  Enumerated  Subsets  and 
Stepwise  Procedures 


Evaluation 

Method 


k= 10 ,  m= 100 


Case 

II 


III 


k=20 ,  m=50 


Case 

II 


III 


Enum.  Subsets 

5 

7 

4 

0 

2 

0 

Stepwise 

3 

______ 

3 

4 

0 

2 

2 

Evaluation 

Method 

=  50,  m: 

Case 

II 

=  20 

III 

k= 

I 

1000, 

Case 

II 

n=  1 

III 

Enum.  Subsets 

”  - - - 

0 

0 

0 

0 

_  _  _ _ 

0 

0 

Stepwise 

it 

II 

H  o 

II 

II 

II 

II 

O  II 

II 

II 

0 

II 

II 

O  II 

II 

II 

II 

II 

O  II 

II 

II 

0 

From  review  of  the  results,  it  was  obvious  that  as  the  number  of 
replications  per  meta-experiment  increased,  both  evaluation  methods 
selected  the  same  subsets.  This  tends  to  validate  the  results  found 
earlier  by  reviewing  the  coverage  and  volume  reduction  figures.  It  is 
reasonable  to  find  that  as  the  coverage  and  volume  reduction  figures 
converge,  the  greater  the  similarity  between  the  variable  subsets 
selected . 


Results  From  Untested  Simulation  Model  Data 

The  results  derived  from  evaluating  the  data/output  from  the 
untested  simulation  model  are  presented  in  Table  4.6.  From  comparing 
the  variable  subsets  selected  by  VSSP  and  SAS,  several  items  were  noted. 
The  most  significant  observation  was  that  regardless  of  selection 
criterion  used,  almost  all  variable  subsets  were  identical.  It  should 


45 


Table  4.6  Effects  of  various  selection  software  and  criterion  on 
selection  of  variable  subsets 

Software 

Package 

Evaluation 

Method 

Selection 

Criterion 

Criterion 

Value 

Variable  Subset* 
Selected 

ft  ft 

UQCP 

Enumerated 

Subsets 

Eq  (62) 

3.156511 

W1  W5  W6  R1  R4  R5  R6 

R7 

v  oor 

Stepwise 

Eq  (62) 

3.156511 

W1  W5  W6  R1  R4  R5  R6 

R7 

*P 

0.306733 

W1  W5  W6  R1  R4  R5  R6 

R7 

Enumerated 

Subsets 

*a 

0.301137 

W1  W5  W6  R1  R4  R5  R6 

R7 

CP 

9.011377 

W3  W4  W5  W6  R1  R4  R5 
R6  R7 

SP 

1 . 196549 

W1  W5  W6  R1  R4  R5  R6 
R7 

Stepwise 

h2 

«p 

0.306733 

W1  W5  W6  R1  R4  R5  R6 

R7 

SAS 

CP 

7.062709 

W1  W5  W6  R1  R4  R5  R6 
R7 

Stepwise 

(Forward 

Selection) 

4 

0.306733 

W1  W5  W6  R1  R4  R5  R6 
R7 

CP 

7.062709 

W1  W5  we  R1  R4  R5  R6 
R7 

Stepwise 

(Backward 

Selection) 

R2 

«p 

0.306769 

W1  W3  W4  W5  R1  R4  R5 
R6  R7 

CP 

9.011541 

W1  W3  W4  W5  R1  B4  R5 
R6  R7 

Stepwise 

(MAXR) 

4 

0.306733 

W1  W5  W6  R1  R4  R5  R6 
R7 

CP 

9.011541 

W1  W4  W5  W6  R1  R4  R5 
R6  R7 

II 

II 

II 

II 

II 

II 

II 

II 

It 

II 

II 

f< 

II 

II 

H 

II 

II 

II 

II 

II 

II 

II 

II 

II 

H 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

H 

«  =  W(i)  is  Work  variable  i,  and  R(i)  is  Routing  variable  i 
•*  =  Using  a  90%  level  of  significance  (i.e.  alpha  «  0.10) 
»••  =  Using  the  standard  SAS  default  significance  criterions 
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also  be  noted  that  of  all  the  other  selection  criterion  used  in  this 

comparison,  Sp  is  the  most  closely  related  to  BCp  and  selected  the  same 

2 

variable  subset.  Also,  the  Ra  criterion  selected  the  same  variable 
subset . 

Other  selection  criterion  used  also  provided  comparable  results. 

2 

The  Rp  criterion  selected  the  same  or  nearly  same  variable  subset 
depending  on  the  point  at  which  further  improvement  was  discarded.  If 
improvement  required  to  be  significant  was  set  to  >  0.001  then  the  same 
variable  subset  was  selected.  However,  if  it  was  set  to  >  0.01  then  a 
subset  without  variable  Rl,  was  selected.  As  noted  in  the  literature 
review,  Chapter  II,  the  level  of  significant  improvement  is  highly 
subjective. 

For  the  Cp  criterion,  only  with  the  enumerated  subsets  procedure 
was  the  best  variable  subset  found.  In  other  procedures,  where  this 
selection  criterion  was  available,  a  near-best  variable  subset  was 
selected.  Regardless  of  these  differences,  each  variable  subset 
selected  contained  one  more  variable  than  that  selected  using  BCp,  and 
differed  only  slightly. 

There  are  two  final  comments  on  the  SAS  derived  results.  First, 
not  all  selection  criterion  were  available  for  all  of  the  SAS  procedures 
employed.  And  second,  the  variable  subset  selected  by  any  procedure,  may 
depend  on  the  level  of  significance  used.  The  default  values  for  the 
SAS  procedures,  used  in  this  evaluation,  were  not  all  set  to  the  same 
level  as  used  in  the  VSSP.  This  may  account  for  the  minor  differences 
noted . 

The  conclusions  drawn  from  these  results  are  summarized  in  the 
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following  chapter.  This  chapter  also  outlines  any  recommendations  for 
further  study  and  research  in  this  area. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

From  the  work  entailed  in  this  thesis,  the  following  conclusions 
and  recommendations  were  derived. 

Conclusions 

From  review  and  evaluation  of  the  results,  the  following 
conclusions  were  drawn  in  regards  to  the  Variable  Subset  Selection 
Program  and  the  performance  of  the  BCp  against  other  common  criterion. 
The  selection  procedures  and  criterion  incorporated  into  the  Variable 
Subset  Selection  Program  performed  as  well  as  expected.  The  results 
obtained  by  either  selection  procedure,  enumerated  subsets  or  stepwise 
(forward  selection),  are  comparable  and  contain  minimal  variations,  even 
when  applied  to  a  small  number  of  replications.  In  addition,  any 
differences  between  the  results  derived  by  either  selection  procedure 
decreases,  as  the  number  of  replications  increase,  until  there  is  no 
difference.  It  was  also  concluded  that  the  known  covariances  between 
the  controls  only  has  a  significant  effect  on  the  coverage  and  volume 
reduction  figures  when  applied  to  a  small  number  of  replications. 
However,  this  advantage  rapidly  disappears  as  the  number  of  replications 
increase . 

In  regards  to  the  performance  of  the  BCp  criterion  in  comparison  to 
other  criterion  in  common  use  today;  the  BCp  criterion  provides 
reasonable  and  comparable  results.  Thus  the  Variable  Subset  Selection 
Program  appears  to  be  a  reasonable  substitute  for  use  by  organizations 
requiring  limited  evaluation  of  this  sort,  where  use  of  a  commercial 
package  may  be  infeasible.  This  infeasibility  may  take  several  forms, 
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notably  expense  or  cost-effect iveness,  limited  access  to  existing 
software,  or  time  involved  in  learning  to  use  the  software  effectively. 

Recommendations 

The  recommendations  derived  from  this  study  fall  into  two 
categories;  1)  further  improvements  to  the  VSSP,  and  2)  further 
study/experimentation  involving  the  BCp  criterion.  There  are  numerous 
improvements  which  can  be  made  to  the  VSSP.  First  and  foremost  is  to 
implement  a  more  efficient  stepwise  (forward  selection)  procedure. 
Primarily  this  entails  developing  a  satisfactory  scheme  for  saving  prior 
pivots  of  the  ‘A‘  matrix  so  subsequent  pivots  are  performed  on  the 
correct  matrix.  The  current  implementation  performs  a  true  stepwise 
(forward  selection)  search,  but  reverts  to  performing  all  pivots  on  each 
subset  evaluated.  This  is  inefficient  and  slows  down  program  execution. 
The  main  benefit  to  be  derived  from  this  recommendation  would  be 
increased  speed  of  evaluation. 

Next,  would  be  the  implementation  of  additional  selection 
procedures,  notably  a  true  stepwise  procedure  would  be  desirable.  Other 
selection  procedures  as  outlined  in  the  Literature  Review  (Chapter  2) 
are  also  viable  candidates. 

Another  recommendation  is  to  revise  the  VSSP  into  a  more  modular 
and  efficient  form.  This  can  be  accomplished  by  breaking  up  the  main 
program  into  subroutines.  Each  selection  procedure  implemented  could  be 
made  a  separate  subroutine  called  by  the  main  program  when  needed. 

Also,  eliminate  duplicate  variables  where  possible  without  affecting 
program  execution.  And  finally,  arrays  and  matrices  shared  by  other 
program  subroutines  could  be  incorporated  into  common  blocks.  This  will 
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eliminate  creation  of  duplicate  constructs  and  cut  down  on  the  amount  of 
memory  required  to  run  the  program. 

Additionally,  there  are  several  recommendations  regarding  further 

study,  analysis,  and  research  using  the  VSSP  and  the  BCp  criterion.  One 

recommendation  is  to  evaluate  output  of  a  simulation  model,  using  the 

VSSP,  to  obtain  the  single  best  variable  subset.  This  could  be  done 

using  both  the  enumerated  and  stepwise  procedures  to  check  for 

consistency  in  the  results.  Next,  run  two  experimental  designs,  one 

using  the  variables  selected  by  the  VSSP  and  the  second  using  all 

variables  in  the  original  full  set.  Then  evaluate  the  experimental 

design  results  using  similar  selection  procedures  but  based  on  other 

2  2 

selection  criterion  (i.e.  Cp ,  B  ,  Ba,  etc.).  Compare  the  final  results 
and  evaluate  the  differences,  if  any.  Does  the  variable  subset  selected 
by  the  VSSP  provide  a  good  starting  point  for  an  experimental  design? 

Another  recommendation  is  to  evaluate  a  set  of  data  while 
increasing  the  number  of  replications,  but  holding  the  number  of  metc.- 
experiments  steady.  How  does  this  affect  the  overall  coverage  and 
volume  reduction  figures?  Can  any  consistent  trends  be  identified?  The 
VSSP  does  not  require  the  use  of  all  data  in  a  datafile  to  be  used  in 
execution. 

And  finally,  experiment  with  data/output  with  known  characteristics 
(i.e.  means  of  controls  and  responses,  and  covariances  between  controls 
and  responses),  to  find  the  point  where  the  results  converge  across  all 
combinations  of  selection  procedures  and  whether  the  covariances  are 
estimated  or  known.  Then  determine  if  this  point  can  be  arrived  at 
analytically.  If  so,  this  type  of  knowledge  could  prevent  gathering 
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more  data  than  necessary  for  an  optimal  evaluation.  This  could  also 
help  contain  costs  associated  with  gathering  data.  Additionally,  this 
may  provide  a  means  of  assessing  reliability  of  results,  depending  on 
the  amount  of  data  available  for  the  evaluation. 
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APPENDIX  A:  FORTRAN  Listing  of  Variable  Subset  Selection  Program 


C 

C  * 
C  * 
C  * 
C  » 

c  « 
c  » 
c  « 
c  * 
c  * 
c  * 
c  » 
c  * 
c  * 
c  » 
c  * 
c  » 
c  « 
c  » 
c  « 
c  « 
c  * 
c  » 
c  * 
c  « 
c  * 
c  « 
c  « 
c  * 
c  * 
c  « 
c  « 
c  » 
c  * 
c  » 
c  * 
c  * 
c  « 
c  » 
c  » 
c  » 
c  * 
c  * 
c  * 
c  » 
c  « 
c  » 
c  » 


PROGRAM  SELECT (INPUT .OUTPUT ,TAPE7 ,TAPE5= INPUT ,TAPE6=0UTPUT) 
*###«*»»#**«*#****#««*#***#**«*«#***«**»*»»**»»#****#«*«»**«#*#»»» 


This  program  provides  both  an  ‘all  possible  regressions* 

(i.e.  enumerated  subsets)  and  a  stepwise  (forward  selection) 
approach  to  selecting  the  best  subset  of  controls  from  a 
given  candidate  set.  It  assumes  that  a  certain  number  of 
meta-experiments  have  been  performed,  each  with  the  same 
number  of  replications.  Once  the  optimal  subset  has  been 
identified,  a  confidence  region  is  constructed  about  the  mean 
vector  for  the  responses.  The  corresponding  coverage  and 
volume  reduction  are  then  tallyed  and  subsequently  summarized. 

This  program  can  be  run  in  two  modes.  The  user  can  either 
estimate  the  covariance  matrix  of  controls  or  incorporate 
it  directly.  The  program  variable  ‘iknow*  dictates  which 
option  is  in  effect  (see  code  below). 

The  program  can  also  be  run  in  the  ‘best  m‘  regressions  mode. 

(  Currently  only  configured  for  estimated  covariance  matrix 
of  controls) 

In  other  words  it  will  compute  the  best  m  subsets  of  each 
possible  subset  size.  This  can  be  of  interest  if  a  single  set 
of  data  is  used. 


PROGRAM  PARAMETERS: 

Z1  =  Max  •  of  candidate  controls  allowed 

Z2  =  Max  *  of  responses  allowed 

Z3  =  Z1  +  Z2 

Z4  =  Max  •  of  best  regressions  which  may  be  kept 

(m  in  ‘m  best*  as  above) 

Z5  =  2»*Z1 

Z6  =  Max  •  replications  per  meta  experiment  allowed 

Z7  =  Max  •  of  meta  experiments  allowed 

Z8  =  (Z3* (Z3+ 1 ) ) / 2 

MAX  =  Maximum  number  of  storage  locations  in  array  A 
for  matrices  created  by  Subroutine  GAUSS. 

CORRESPONDING  PRIMARY  VARIABLES,  INITIALIZED  BY  USER: 

NX  *  t  of  candidate  controls 

NY  =  •  of  responses 

NVAR  =  NX  ♦  NY 

KEEPERS  =  •  of  best  regressions  to  be  kept 

(m  in  'm  best'  as  above) 
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o  n  o  o  o  o  o  o  an  o  o  o  o  o  o 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  SELECT 

•««««**»«««««*»««*«*«*»«««*«*»»«»«•««»««»«**««*««**•*«»«««•««* 
*  PROGRAM  PARAMETERS  AND  VARIABLE  INITIALIZATION  • 

»««««««*««»««««*•**««»•«««•««**•«•*«»«««»•»••«««•»*•*««»**#**« 

INTEGER  Zl,  Z2 ,  Z3,  Z4 ,  Z5 ,  Z6 ,  Z7 ,  Z8 ,  MAX 

PARAMETER  (Zl=8,  Z2=8,  Z3=Z1+Z2,  Z4=6,  Z5=2*«Z1) 

PARAMETER  (Z6=50,  Z7=100,  Z8= ( (Z3» (Z3+ 1 ) 5 /2) ) 

PARAMETER  (MAX=50) 

COMMON  /BLK1/  SIG,  KK,  IQQ,  IP 

CHARACTER  TITLE*25,  RESPONS (Z2) *25 ,  CONTROL (Zl) *25 
CHARACTER  INFILE»25,  OUTFILE*25,  COVFILE*25 
CHARACTER  ANSWER,  XFILE»25 

INTEGER  NK(Z3),  MODELS (Z5 ,Z1) ,  IBUFF(Zl),  IX,  KK,  IP 
INTEGER  IH(Z3),  ICOVER(4) ,  ICTOT(4),  NBR(6) ,  IIN,  IQQ 
INTEGER  NX,  NY.  NVAR,  KEEPERS,  KNX,  NUMBERS,  META 
INTEGER  I KNOW,  IWRITE,  METHOD,  II,  12,  IZ,  KP,  K 
INTEGER  TMV,  INDKZl  +  l),  IND2  (Z1  + 1 )  ,  TIND  (Zl  + 1 ) 

REAL  COVCV(Zl ,Z1) ,  VECMUC(Zl),  CBAR(Zl),  VECCBAR(Zl) 

REAL  FF(0:Z1),  WKAREA(Z2) ,  RSS(Z2,Z2),  DUM(Z2) 

REAL  TARGET (Z2.Z2) ,  VECYBAR(Z2) ,  VECMUY(Z2) ,  YBARCZ2) 

REAL  A(MAX,Z3,Z3) ,  VCV(Z8),  FULL(Z3,Z3) 

REAL  REGR(Z4 ,Z1 ,2) ,  BUFF(Z4) ,  BUFF2(Z4) 

REAL  VR(2)  ,VOLRED(2)  , COVERAGE)  ,  X(Z6,Z3),  TEMP(Z3) 

REAL  XM(Z3) ,  SUMDEV(2) ,  SUMVU(2) ,  SiG 


KNX  =  2*»NX 

NUMREPS  =  *  replications  per  neta  experiment 
META  =  i  of  meta  experiments 

NOTE: 

IN  SUBROUTINE  COVER  :  Zl  AND  Z2  IN  THE  PARAMETER 

STATEMENT.  MUST  BE  SET  TO 
Zl  AND  Z2,  RESPECTIVELY,  OF 
THE  MAIN  PROGRAM  PARAMETER 
STATEMENT 


»  BRIEF  PROGRAM  INTRODUCTION  AND  INITIAL  DATA  INPUT/OUTPUT  « 
•  ROUTINE  • 
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o  o  o  o  o 


c  ******************************************* 

C  *  PROGRAM  INTRODUCTION  * 

C  ******************************************* 

c 

PRINT  « ,  ‘THIS  IS  THE  VARIABLE  SUBSET  SELECTION  PROGRAM' 
PRINT  »,  •' 

PRINT  *,  '  THIS  PROGRAM  IS  DESIGNED  TO  EVALUATE  OUTPUT' 
PRINT  »,  'FROM  A  SIMULATION  MODEL  AND  DETERMINE  THE  OPTIMAL 
PRINT  *,  'SUBSET  OF  VARIABLES,  BASED  ON  THE  "BEST  CONTROLS 
PRINT  *,  'CRITERION.  YOU  MAY  CHOOSE  EITHER  AN  ENUMERATED' 
PRINT  »,  ‘SUBSETS  OR  STEPWISE  (FORWARD  SELECTION)  APPROACH. 
PRINT  »,  " 

PRINT  *,  " 

utiitotmmtmoitiitoHiimiHHHi 

*  DATA  INPUT  » 

******************************************* 

10  PRINT  *,  ‘DO  YOU  WISH  TO  ENTER  PROGRAM  DATA  MANUALLY  OR  BY' 
PRINT  *,  'DATAFILE?  ENTER  M  OR  D:' 

READ  »,  ANSWER 

IF  ((ANSWER. EQ. ’ D ' ) .OR. (ANSWER. EQ. *d’))  THEN 
PRINT  »,  'ENTER  NAME  OF  THE  INPUT  DATAFILE:  ' 

READ  *,  INFILE 

OPEN (UNIT= 10 ,  FILE=INFILE ,  STATUS= ’ OLD 1 ) 

READ (10,*)  NX.  NY,  KEEPERS,  NUMREPS ,  META 
READ (10,*)  I KNOW,  IWRITE,  METHOD 
READ (10,*)  SIG 

READ (10,*)  (VECMUC(I) ,1=1, NX) 

READ(10,»)  (VECCBAR(I) ,1=1 ,NX) 

READ(10,»)  (VECMUY(I) , 1=1 ,NY) 

READ(10,»)  (VECYBAR(I) ,1=1 ,NY) 

READ (10,*)  TITLE 
READdO,*)  (CONTROL (I)  ,1  =  1, NX) 

READdO,*)  (RESPONS(I)  ,1=1  ,NY) 

IF  (IKNOW.EQ. 1)  THEN 
READdO,*)  COVFILE 
IF  (COVFILE. NE. INFILE)  THEN 

OPEN (UNIT5 15 ,  FILE=COVFILE ,  STATUS5 ’ OLD’ ) 

DO  15  1=1, NX 

READ (15,*)  (COVCV(I,J) ,J=1,NX) 

15  CONTINUE 

CLOSE (15) 

ELSE 

DO  20  1=1, NX 

READdO,*)  (COVCV(I.J)  ,J=1,NX) 

20  CONTINUE 

END  IF 
END  IF 

READdO. »)  XFILE 

OPEN (UNIT=20 ,  FILE=XFILE ,  STATUS5 ’OLD’ ) 

ELSE 
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25 


30 


35 


IF  ( (ANSWER. NE.'M'). AND. (ANSWER. NE.’m’))  THEN 
PRINT  *,  "INVALID  INPUT.  TRY  AGAIN." 

PRINT  *.  "" 

GOTO  10 


ENDIF 

PRINT  » 

PRINT  * 

PRINT 

PRINT 

READ  « 

PRINT 

READ  * 

PRINT 

PRINT 

READ  * 

PRINT 

PRINT 

READ  * 

PRINT 

PRINT 

READ  * 

PRINT 

PRINT 

READ  » 

PRINT 

PRINT 

READ  * 

PRINT 

PRINT 

READ  * 

PRINT 

PRINT 

READ  * 

PRINT 

PRINT 


"ENTER  THE  FOLLOWING  VARIABLE  VALUES:" 

•  • 

"  INPUT  •  OF  CANDIDATE  CONTROLS  (PROGRAM  * 

"  LIMIT  =  " ,21 , " ) : * 

NX 

"  INPUT  •  OF  RESPONSES  (PROGRAM  LIMIT  =  \Z2,"): 
NY 

*  INPUT  •  OF  BEST  REGRESSIONS  TO  KEEP" 

"  (PROGRAM  LIMIT  =  ",Z4,"):  * 

KEEPERS 

*  INPUT  *  OF  REPLICATIONS  PER  META  EXPERIMENT" 

*  (PROGRAM  LIMIT  =  " ,Z6,") :  * 

NUMREPS 

"  INPUT  «  OF  META  EXPERIMENTS  DESIRED  * 

"  (PROGRAM  LIMIT  =  " ,Z7, ") :  * 

META 

"  INPUT  WHETHER  COVARIANCE  MATRIX  OF  CONTROLS  IS" 
"  ESTIMATED  (0),  OR  KNOWN  (1):  " 

I  KNOW 

*  INPUT  WHETHER  YOU  WANT  THE  META  EXPERIMENT  MODE 

"  (0).  OR  THE  BEST  M  REGRESSIONS  MODE  (1):  * 

IWRITE 

"INPUT  WHETHER  YOU  WANT  TO  USE  ENUMERATED  SUBSETS* 
"(0),  OR  STEPWISE  [FORWARD  SELECTION]  (1)  METHOD: 
METHOD 

"  INPUT  LEVEL  OF  SIGNIFICANCE  OF  TEST  (e.g.  " 

"  90X  =  0.90) :  " 

SIG 

■  IN?;T  THE  KNOWN  MEAN  FOR  EACH  CONTROL:  " 


DO  25  1=1. NX 

PRINT  «.  "  VECMUCC  ,1,")  =  " 

READ  »,  VECMUC(I) 

CONTINUE 

PRINT  «.  "  INPUT  AVERAGE  OF  INPUTS  FOR  EACH  CONTROL:  " 
PRINT  »,  "" 

DO  30  1=1, NX 

PRINT  »,  "  VECCBAR ( ’ , I , ’ )  =  " 

READ  «,  VECCBAR ( I ) 

CONTINUE 

PRINT  •,  "  INPUT  ESTIMATED  MEAN  FOR  EACH  RESPONSE:  " 
PRINT  »,  *" 

DO  35  1=1, NY 

PRINT  «,  "  VECMUYC  ,1.")  =  " 

READ  *.  VECMUY ( I ) 

CONTINUE 
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PRINT  «,  •  INPUT  AVERAGE  OF  CONTROLLED  OBSERVATIONS  FOR 
PRINT  «,  '  EACH  RESPONSE:  ' 

PRINT  «,  '• 

DO  40  1=1, NY 

PRINT  *,  -  VECYBARC  ,1 , ')  =  * 

READ  *,  VECYBAR(I) 

40  CONTINUE 

PRINT  *,  '  INPUT  ANALYSIS  TITLE:  * 

READ  » ,  TITLE 

PRINT  »,  ‘  INPUT  NAMES  OF  CANDIDATE  CONTROLS:  ' 

PRINT  *,  ** 

DO  45  1=1, NX 

PRINT  *,  '  CONTROLC  ,1,')  =  ' 

READ  »,  CONTROL (I) 

45  CONTINUE 

PRINT  «,  ‘  INPUT  NAMES  OF  RESPONSES:  ’ 

PRINT  *,  " 

DO  50  1=1, NY 

PRINT  «,  '  RESPONSE  C,  I ,’ )  =  * 

READ  *,  RESPONS ( I ) 

50  CONTINUE 

IF  (IKNOW.EQ.l)  THEN 

55  PRINT  »,  ‘WILL  YOU  ENTER  KNOWN  COVARIANCE  MATRIX,  OF 

PRINT  *.  'CONTROLS,  MANUALLY  OR  BY  DATAFILE  (M  or  D) : 
READ  #,  ANSWER 

IF  ((ANSWER. EQ.'D'). OR. (ANSWER. EQ.'d'))  THEN 

PRINT  *,  'ENTER  NAME  OF  DATAFILE  CONTAINING  KNOWN 
PRINT  ».  'COVARIANCE  MATRIX  OF  CONTROLS:  ' 

READ  «,  COVFILE 

OPEN(UNIT= 15 ,  FILE=COVFILE ,  STATUS= ’OLD ’ ) 

DO  60  1=1, NX 

READ ( 15 , »)  (COVCV(I ,J) ,J=1 ,NX) 

60  CONTINUE 

CLOSE (15) 

ELSE 

IF  ( (ANSWER. NE. 'M') .AND. (ANSWER. NE. 'm'))  THEN 
PRINT  «,  'INVALID  RESPONSE,  TRY  AGAIN.' 

GOTO  55 
END  IF 

PRINT  *,  '  INPUT  THE  REQUESTED  VALUES:  ' 

PRINT  *,  " 

DO  65  1=1, NX 
DO  65  J= 1 ,NX 

PRINT  *,  '  COVCVC.I,  V.J,*)  =  ' 

READ  »,  COVCV(I.J) 

65  CONTINUE 

END  IF 
END  IF 

70  PRINT  *,  'WILL  YOU  ENTER  THE  [CONTROLS  1  RESPONSES ]  DATA  ' 
PRINT  «,  'MANUALLY  OR  BY  DATAFILE  (M  or  D) :  ' 

READ  »,  ANSWER 
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noon 


IF  ( (ANSWER. EQ. ’D‘ ) .OS. (ANSWER. EQ. ‘d‘ ) )  THEN 

PRINT  «,  ‘ENTER  NAME  OF  [CONTROLS ! RESPONSES ]  DATAFILE:  ‘ 
READ  *,  XFILE 

OPEN (UNIT=20 ,  FILE=XFILE ,  STATUS5 'OLD' ) 

ELSE 

IF  (( ANSWER. NE. ‘M‘) .AND. (ANSWER. NE. ‘m‘) )  THEN 
PRINT  »,  ‘INVALID  RESPONSE,  TRY  AGAIN.’ 

GOTO  70 
ENDIF 

PRINT  *,  ‘  MATRIX  VALUES  WILL  BE  REQUESTED  AS  REQUIRED’ 
ENDIF 

75  PRINT  «,  ‘DO  YOU  WANT  YOUR  INPUTS  SENT  TO  A  DATAFILE  FOR  * 
PRINT  »,  ‘FUTURE  USE  (Y/N) :  ‘ 

READ  *,  ANSWER 

IF  ((ANSWER. EQ.‘Y‘) .OR. (ANSWER. EQ.‘y‘))  THEN 

PRINT  *,  ‘  ENTER  NAME  OF  INPUT  DATAFILE  TO  CREATE:  * 

READ  «,  OUTFILE 

OPEN (UNIT=25 ,  FILE=OUTFILE) 

WRITE ( 25 ,« )  NX , NY , KEEPERS , NUMREPS , META 
WRITE (25 , *)  IKNOW, I WRITE .METHOD 
WRITE (25 , *)  SIG 
WRITE (25 , »)  (VECMUC(I) ,1=1 ,NX) 

WRITE (25 , *)  (VECCBAR(I) ,1=1, NX) 

WRITE (25 , *)  (VECMUY(I) ,1=1 ,NY) 

WRITE (25 , » )  (VECYBAR(I) ,1=1 ,NY) 

WRITE (25 ,  *)  TITLE 

WRITE (25 , *)  (CONTROL ( I ) ,1=1, NX) 

WRITE (25 , *)  (RESPONS ( I ) ,1=1 ,NY) 

IF  ( IKNOW. EQ.l)  THEN 
WRITE (25 , »)  COVFILE 
IF  (COVFILE. EQ. OUTFILE)  THEN 
DO  80  1=1, NX 

WRITE (25 , »)  (COVCV(I ,J) ,J=1 ,NX) 

80  CONTINUE 

ENDIF 
ENDIF 

WRITE (25, »)  XFILE 
CLOSE (25) 

ELSE 

IF  ( (ANSWER. NE. ‘N‘) .AND. (ANSWER. NE.‘n‘))  THEN 
PRINT  »,  ‘  INVALID  INPUT,  TRY  AGAIN.’ 

GOTO  75 
ENDIF 
ENDIF 
ENDIF 
C 

NVAR  =  NX  ♦  NY 
KNX  =  2*«NX 


»  TEST  PRIMARY  VARIABLES  INPUT  » 
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c 

1  =  0 

PRINT  »,  " 

IF  (NX.GT.Z1)  THEN 
PRINT  1575,  Z1 
1  =  1  +  1 
ENDIF 

IF  (NY.GT.Z2)  THEN 
PRINT  1580,  Z2 
1  =  1  +  1 
ENDIF 

IF  ( KEEPERS. GT.Z4)  THEN 
PRINT  1585 
PRINT  1586,  Z4 
1  =  1  +  1 
ENDIF 

IF  (NUMREPS . GT . Z6)  THEN 
PRINT  1590 
PRINT  1591,  Z6 
1  =  1  +  1 
ENDIF 

IF  (META.GT.Z7)  THEN 
PRINT  1595,  Z7 
1  =  1  +  1 
ENDIF 

IF  (I.GT.O)  THEN 
PRINT  1600,  I 
STOP 
ENDIF 
C 

C  a*###*#*##*###***##*##*##**#*#*###***#####* 

C  *  INITIAL  DATA  OUTPUT  » 

C  »*»#******»»*»»###**»***»»******#»»*******» 

c 

PRINT  »,  ’ENTER  NAME  OF  FILE  FOR  PROGRAM  OUTPUT:  * 
READ  »,  OUTFILE 
PRINT  «,  ’’ 

OPEN (UNIT=30 ,  FILE=OUTFILE ,  STATUS= ’ NEW’ ) 

WRITE(30, 1500)  TITLE , META , NUMREPS , META»NUMREPS 
WRITEOO,  1515) 

WRITE(30, 1505)  META»NUMREPS 
DO  90  1=1, NY 

WRITE (30, 1510)  I , RESPONS ( I ) .VECYBAR(I) ,VECMUY(I) 
90  CONTINUE 

WRITE(30, 1515) 

WRITE(30, 1520)  META*NUMREPS 
DO  95  1=1, NX 

WRITEOO,  1510)  I , CONTROL(I)  ,VECCBAR(I)  ,VECMUC(I) 
95  CONTINUE 

WRITEOO,  1515) 

C 
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C  «  DECLARE  IF  COVARIANCE  MATRIX  « 

C  »  OF  CONTROLS  IS  KNOWN  OR  ESTIMATED  » 

C  «  « 

C  *  IKNOW  *  0,  COV  MATRIX  ESTIMATED  » 

C  *  IKNOW  =  1,  KNOWN  COV  MATRIX  USED  » 

C  »**#»##*»#**»##***»*t»«#t***»#*»#*#**» 

C 

IF  ( IKNOW. EQ.O)  THEN 
WRITE (30, 1525) 

ELSE 

WRITE (30 , 1530) 

END  IF 

WRITE(30, 1515) 

C 

C  tiiKUHimiMmitmHotiitdiitx* 

C  »  PROVIDE  HEADING  FOR  REMAINDER  OF  * 

C  «  PROGRAM  OUTPUT  « 

C  »#»#»*»«#***«*»»*»**»»§»»§##*****»#»#» 

c 

WRITE(30, 1515) 

WBITE(30, 1535) 

WRITE (30 ,1515) 

C 

C  tHHIMIIHIXiMmOHIimitl 

C  «  DEFINE  INPUT  VECTOR  FOR  IMSL  * 

C  *  SUBROUTINE  'BECOVM'  » 

C  »*»»»#»*«»»*»»»*»*»***»»#*»*»#»** 

c 

NBR ( 1 ) =Z3 
NBR(2) =NUMREPS 
NBR (3) =NUMREPS 
NBR (4) *  1 
NBR (5) = 1 
NBR (6) =0 
IX=Z6 
C 

C  •  BEGIN  MAIN  PROGRAM  • 

C  . . . 

C 

C  *  MAKE  THE  F  TABLE  » 

C 

IP=NY 

KK=NUMREPS 

CALL  FTABL(FF,NX,Z1) 

PRINT  «,’F  TABLE: ’ , (FF ( I ) ,1=0, NX) 

PRINT  *,  ” 

C 
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o  o  o  o  o  o  o  o  o  o 


C  MHmmtoMommtmmmiHiiim* 

C  *  INITIALIZE  COVERAGE  AND  VOLUME  « 

C  »  REDUCTION  ACCUMULATORS  * 

C  »«»#»»**#**#»»**»»»»»»*»*#»*»«##«»*#»»**#»« 

c 

DO  100  IZ=  1 , 4 
ICTOT ( IZ) =0 
100  CONTINUE 
C 

DO  105  IZ= 1 . 2 
SUMDEV ( IZ) =0 . 

SUMVU ( IZ) =0 . 

VR(IZ) =0. 

105  CONTINUE 

*•«*«**•«•««*««««««*«**«****•*«*«•«**««•«»*«»*»• 

•  THIS  IS  THE  META  EXPERIMENT  LOOP  « 

**#*###***###*#***«#*«##***#**«**«##*#*###**#**# 

DO  1000  MM= 1 , META 

##»»##*»»»#«**#**»#*»*#*###*»#*#« 

*  INITIALIZE  ARRAYS  « 

*«*«*»«**«*«*««««f«**«»**««««««*« 

NUMREG=0 

DO  110  IZ- 1 .KEEPERS 
DO  110  JZ= 1 ,NX 
DO  110  KZ= 1 , 2 

REGR(IZ, JZ.KZ) =0. 

110  CONTINUE 

C 

DO  115  IZ= 1 , KNX 
DO  115  JZ= 1 ,NX 
MODELS (IZ. JZ) =0 
115  CONTINUE 

C 

DO  120  IZ= 1 ,NVAR 
DO  120  JZ= 1 ,NVAR 
DO  120  KZ= 1 ,NVAR 
AdZ,  JZ.KZ) =0. 

120  CONTINUE 

C 

DO  125  IZ=1, KEEPERS 
BUFF (IZ) =0 
BUFF2 (IZ) =0 
125  CONTINUE 

C 

DO  130  IZ=1,NX 
IBUFF ( IZ) *0 
130  CONTINUE 


C 


c  *»*»*»»»*»»**###»»»«#«*#*•******»»*#*#*»#«» 


C  •  BEAD  THE  DATA  * 
C  *  (EACH  RECORD  =>  [CONTROLS ! RESPONSES ] )  » 
C  *  COMPUTE  THE  COVARIANCE  MATRIX  « 
C  *  SAVE  SAMPLE  MEANS  » 
C  *  BOUND  THE  GENERALIZED  VARIANCE  « 


C  #»**»***#**»##*t»#*»*#»»****»#*»*»*»»#**»»# 

c 

IF  ((ANSWER. EQ. ’M') .OR. (ANSWER. EQ. ’m’))  THEN 

PRINT  »,  'ENTER  ELEMENTS  OF  [CONTROL! RESPONSE]  MATRIX' 
PRINT  *,  '' 

DO  131  1=1 .NUMREPS 
DO  132  J= 1 , NVAR 

PRINT  «,  '  XC.I.'.'.J,')  =  ' 

READ  »,  X(I,J) 

132  CONTINUE 

131  CONTINUE 

ELSE 

DO  135  1=1, NUMREPS 

READ (20, «) (X(I,J) , J=1 ,NVAR) 

135  CONTINUE 

ENDIF 
C 

CALL  BECO VM ( X , I X , NBR , TEMP , XM , VC  V , I ER ) 

C 

DO  140  1=1, NX 
CBAR(I) =XM(I) 

140  CONTINUE 

C 

DO  145  I=NX+1,NVAR 
YBAR ( I -NX) =XM ( I ) 

145  CONTINUE 

C 

CALL  VCVTSF (VCV ,Z3 .FULL ,Z3) 

C 

DO  150  1=1, NVAR 
DO  150  J= 1 , NVAR 

A ( 1 , I , J ) =FULL(I , J) 

150  CONTINUE 

C 

IS*  1 

DO  155  11=1, NY 
DO  155  JJ=1 ,NY 

IF  (JJ.GE.II)  THEN 

RSS(II ,JJ)=A(IS,NX+II ,NX+JJ) 

RSS (JJ , II) =RSS (II , J  J ) 

ENDIF 

155  CONTINUE 

C 

CALL  LINV3F (RSS ,DUM,4 , NY,Z2 ,D1 ,D2 , WKAREA , IER) 

IF  (IER.NE.O)  PRINT  *,'I  DIED  BELOW  155  (MAIN)' 
DET=D1*2**D2 
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OOOOOCi  O  O  O  O  O  O  O  OOOOOOO  O  CJ>  O  O  O  O  C>  O  O  CJ 


BIQ= (FLOAT (NUMREPS-1) /FLOAT (NUMREPS-NX-2) ) «»NY 
TW0=2»BIG»D1»2«»D2 


«***f«*****t**»«««»**»tt«*«**tt**tt***«»*t* 

«  STUFF  THE  BOOKKEEPING  ABRAY  WITH  THE  » 
«  BOUND  » 


DO  160  11=1 .KEEPERS 
DO  160  JJ=1 ,NX 

REGR ( I I , JJ , 1 ) =TWO 
160  CONTINUE 

*«««****«*««•«*•«•«««««*»»»«»•••»*««*»«**««»*««« 

*  THE  FOLLOWING  SECTION  IS  PERFORMED  IF  THE  « 

»  METHOD  OF  ENUMERATED  SUBSETS  IS  SELECTED  * 
»  (i.e.  METHOD  =0)  » 

IF  (METHOD. EQ.O)  THEN 

«»«•««•«•«*««•*•••*«««•«*«**•««**»«*•»***** 

*  CONDUCT  BINARY  SEARCH  OF  THE  * 

*  REGRESSION  TREE:  « 

«  FURNIVAL  AND  WILSON  (1974)  » 

##*#**»***#*»»»*»*##»***»#*###*»***#*###*#» 


K=NX 


DO  165  L= 1 , K 
NK (L) =0 

165  CONTINUE 


NK (K+ 1 ) = 1 
L=  1 

170  NK (L) * J 


DO  175  M=L,K 

IF  (NK(M+1) .EQ. 1)  GOTO  180 
175  CONTINUE 


180  IB=K-M+ 1 

IS=K-L+2 
IP=K-L+1 
KP=NVAR 

CALL  GAUSS (IB,IS,IP,A,KP, MAX , Z3 ) 


*  CALCULATION  OF  THE  GENERALIZED  RESIDUAL  * 
t  COVARIANCE  « 
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DO  185  11=1, NY 
DO  185  JJ=1,NY 

IF  (JJ.GE.II)  THEN 

RSS (II , JJ) =A(IS , NX+II , NX+JJ) 

RSS (JJ,II)=RSS(II,JJ) 

END  IF 

185  CONTINUE 

C 

IF  ( I KNOW. EQ . 0)  THEN 

CALL  LINV3F (RSS ,DUM,4 ,NY,Z2 ,D1 ,D2 , WKAREA, IER) 

IF  (IER. NE . 0)  PRINT  *,'I  DIED  BELOW  185' 
DET=D1«2»*D2 
END  IF 
C 

C  «*»«#»»*»»##*»##*****#«»****«»«*»»»»»#***** 

C  *  BOOKKEEPING  LOGIC  TO  SAVE  M=KEEPERS  « 

C  «  BEST  REGRESSIONS  OF  ALL  J  SUBSETS  SIZES  » 

C  **#»»**»*»*#***#******o**##***«*****»***»» 

c 

MV=0 

DO  190  N=1 ,NX 
MV=MV+NK  (N) 

190  CONTINUE 

C 

IF  (IKNOW.EQ.O)  THEN 

CONST  = ( FLOAT ( NUMREPS -  1 ) / FLOAT ( NUMREPS - MV - 1 ) ) 
DET=DET*CONST«*NY 
ELSE 

CALL  COVKNOW ( RSS , NY , Z2 , FULL , NVAR , Z3 , TARGET , DUM , 
&  NUMREPS, MV, DET) 

END  IF 
C 

DO  195  J=l, KEEPERS 

IF  (DET.LT.REGR(J ,MV, 1) )  THEN 
NUMREG*NUMREG+ 1 
DO  200  JJ=J, KEEPERS- 1 

BUFF (JJ+ 1 ) =REGR ( JJ , MV , 1 ) 

BUFF2 (JJ+1) =REGR ( J J , MV , 2) 

200  CONTINUE 

REGR(J,MV, 1) =DET 
REGR ( J , MV , 2) =NUMREG 
DO  205  JJ=J+1, KEEPERS 
REGR( JJ , MV , 1 ) =BUFF (JJ) 

REGR ( JJ ,MV , 2) =BUFF2 ( JJ) 

205  CONTINUE 

CALL  KEEPIT( NUMREG , NK , NX , MODELS , Z 1 , Z3 , Z 5 ) 
GOTO  210 
END  IF 

195  CONTINUE 

210  CONTINUE 

C 
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DO  215  L= 1 , K 

IF  (NK(L).EQ.O)  GOTO  170 
NK (L) “0 

215  CONTINUE 

END  IF 
C 

C  HommtmimmmmtittotHtiHtoimiti 
C  *  ENUMERATED  SUBSETS  CODE  ENDS  ON  ABOVE  LINE  » 
C  imtoHomoimtmiotmmmoxomH 
C 

c  #»»**#»#**##**#*»***§#»»*»«****»#»****###*«*»*#* 
C  «  THE  FOLLOWING  SECTION  IS  PERFORMED  IF  THE  « 
C  *  STEPWISE  PROCEDURE  [FORWARD  SELECTION]  IS  * 
C  *  SELECTED  (i.e.  METHOD  =1)  « 

c 


IF  (METHOD. EQ. 1)  THEN 
C 

C  ****#»**#*#**»*#»####*##»»#*#*»#»*«»*#«*»** 


C  »  CONDUCT  STEPWISE  (FORWARD  SELECTION)  * 
C  *  SEARCH  OF  THE  REGRESSION  TREE.  « 
C  *  THIS  IS  A  MODIFIED  VERSION  OF  THE  » 
C  *  NATURAL  SEARCH  PROCEDURE  WRITTEN  BY:  » 
C  «  FURNIVAL  AND  WILSON  (1974)  » 


c 


K=NX 

C 

DO  220  L=1 ,K+1 
IND1 (L) =0 
IND2 (L) =0 
TIND (L) =0 

220  CONTINUE 


M=K 
IB=0 
IS=  1 


TMV=  1 


C 

225  IB=MOD ( IB , MAX) ♦  1 

C 

DO  230  L=M,K 

IF  (IND2(L) . LT.L)  GOTO  230 
IND2 (L-1)=IND2(L-1)+1 
IND2 (L) =IND2 (L- 1 ) 

230  CONTINUE 


C 

235 


C 


IND2 (K) =IND2 (K) + 1 
IS=MOD(IS , MAX) *  1 


MV=0 
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o  o  o  o  o  o 


DO  240  1 1  =  1 tK 

IF  ( IND2  (ID  . GT . 0)  MV=MV*1 
240  CONTINUE 

C 

IF  (MV.GT.TMV)  THEN 
TMV=TMV+ 1 
DO  245  11  =  1, K 

IND1  (ID  =TIND  (1 1 ) 

245  CONTINUE 

END  IF 
C 

IF  (MV.EQ.l)  GOTO  260 
C 

DO  250  11  =  1, K 
FLAG=0 

DO  255  12=1, K+l 

IF  ( IND2  (12)  .EQ.INDKIl))  FLAG=1 
255  CONTINUE 

IF  (FLAG.EQ.O)  GOTO  295 
250  CONTINUE 

C 

260  CONTINUE 
IP=IND2 { K ) 

KP=NVAR 
IB2=  1 
IS2=2 

IF  (MV.GT.2)  THEN 

DO  261  11=1, NX 

IF  (IND2  (ID  .GT.O)  THEN 
IF  (IS2.EQ.2)  THEN 
IS2=3 
ELSE 
IS2=2 
END  IF 

IP=IND2 (1 1) 

CALL  GAUSS (IB2,IS2,IP,A, KP ,MAX ,Z3) 
IB2=IS2 
END  IF 

261  CONTINUE 
ELSE 

CALL  GAUSS (IB, IS, IP, A, KP, MAX ,Z3) 

IS2=IS 
END  IF 


«  CALCULATION  OF  THE  GENERALIZED  RESIDUAL  * 
«  COVARIANCE  » 
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DO  265  11=1, NY 
DO  265  JJ=1,NY 

IF  (JJ.GE.II)  THEN 

RSS (II , JJ) =A(IS2 , NX+II ,NX+JJ) 

RSS (JJ , II) =RSS ( I I , JJ) 

END  IF 

265  CONTINUE 
C 

IF  (IKNOW.EQ.O)  THEN 

CALL  LINV3F (RSS ,DUM,4 , NY.Z2 ,D1 , D2 , WKAREA , IER) 

IF  (IER.NE.O)  PRINT  *,'I  DIED  BELOW  265’ 
DET=D1*2»*D2 
END  IF 
C 

C  »»*#*#»***»»*»t#***»»*»*«*#*t##«**f»*»***## 

C  »  BOOKKEEPING  LOGIC  TO  SAVE  M= KEEPERS  * 

C  »  BEST  REGRESSIONS  OF  ALL  J  SUBSETS  SIZES  » 

C  *t»*#»#*»**#**#t##»##»****#»»*»»»»»***###** 

c 

IF  (IKNOW.EQ.O)  THEN 

CONST= (FLOAT (NUMREPS- 1 ) /FLOAT (NUMREPS- MV- 1 ) ) 
DET=DET»CONST»*NY 
ELSE 

CALL  COVKNOW ( RSS , NY , Z2 , FULL , NVAR , Z3 , TARGET , DUM , 
&  NUMREPS, MV, DET) 

END  IF 
C 

DO  266  11=1, NX 
NK(I 1) =0 

266  CONTINUE 
NK(NX+1)=1 

C 

DO  270  J=l, KEEPERS 

IF  (DET . LT . REGR ( J , MV , 1 ) )  THEN 
NUMREG=NUMREG+ 1 
DO  275  JJ=J, KEEPERS- 1 

BUFF (JJ+ 1 ) =REGR ( JJ , MV , 1 ) 

BUFF2 (JJ+ 1 ) =REGR ( JJ , MV , 2) 

275  CONTINUE 

REGR ( J , MV , 1 ) =DET 
REGR ( J , MV , 2 ) *  NUMREG 
DO  280  JJ=J+1 .KEEPERS 
REGR( JJ ,MV , 1 ) =BUFF ( JJ) 

REGR( JJ ,MV , 2) =BUFF2 (JJ) 

280  CONTINUE 

DO  285  1 1  =  1, NX 

TIND(I 1) =IND2(I 1) 

IF  ( IND2  (ID  .  GT.  0)  THEN 
NK  (NX+ 1 -IND2 ( 1 1 ) )  =  1 
END  IF 

285  CONTINUE 

CALL  KEEP I T ( NUMREG , NK , NX , MODELS , Z 1 , Z3 , Z  5 ) 
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oooooo 


GOTO  290 


ENDIF 

270 

CONTINUE 

290 

CONTINUE 

295 

CONTINUE 

IF  (IMD2 (K) . LT.K)  GOTO  235 
IS=IS-1 

IF  (IND2 (M) . EQ. M)  M=M-1 
IF  (M.GT.O)  GOTO  225 
ENDIF 
C 

C  »»*»#****#**»»**«*t##tt*««»**»#*»»*#»***#**#*#*** 

C  «  STEPWISE  (FORWARD  SELECTION)  CODE  ENDS  ON  • 

C  «  ABOVE  LINE  * 

C  IHmiOiimOOMmiMOtHHHHHIHHHH 
C 

C  »*#**»**»**#»tt»#*»»*»t**»***#**#«*»*#*tt**#» 

C  »  THIS  BLOCK  IS  FOR  BEST  M  SUBSETS  MODE  « 

C  *  OF  OPERATION  « 

C  HtoHmmoHtiHtHomtommmH 

C 

IF  (IWRITE.EQ. 1)  THEN 
DO  300  1=1 ,NX 

WRITE (30 ,1540)  KEEPERS , I 
DO  300  J=1 .KEEPERS 
I VAR=0 
I  IN=0 

DO  305  II=NX, 1,-1 
IVAR=IVAR+ 1 

IF  (IFIX(REGR ( J , I , 2) ♦ . 0001) .EQ . 0)  GOTO  300 
IF  (MODELS (IFIX (REGR (J , I , 2 ) + . 0001) ,11). EQ . 1 )  THEN 
IIN=IIN+1 
IBUFF (IIN) =IVAR 
ENDIF 

305  CONTINUE 

RDET=REGR(J ,1,1) 

WRITE (30, 1545)  MM, RDET , ( IBUFF ( I J) , IJ= 1 . IIN) 

300  CONTINUE 

ENDIF 


»  FOR  EACH  SUBSET  COMPUTE  THE  CRITERION  » 
*  AND  SAVE  THE  MINIMUM  * 


IF  ( I WRITE . EQ . 0)  THEN 
IP=NY 

KK=NUMREPS 
DO  310  IQ= 1 ,NX 
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o  o  o  o  o 


IF  (IKNOW.EQ.O)  THEN 

REGR ( I , IQ , 1 ) =REGR ( 1 , IQ , 1 ) «C3 (KK , IQ , IP) « 

&  CFRONT (KK,IQ,IP)»FF(IQ) *C5 (KK , IQ , IP) 

ELSE 

REflRd.IQ,  l)*REflR(I  , IQ, 1) *C4 (KK , IQ. IP) # 

&  CFRONT (KK,IQ,IP)»FF(IQ) »C5 (KK , IQ » IP) 

ENDIF 

IF  (IQ.EQ.l)  RMIN=REGR ( 1 , IQ , 1 ) *1000 . 

310  CONTINUE 

DO  315  IQ=1,NX 

IF  (REGR (1, IQ, 1) .LT.RMIN)  THEN 
RMIN=REGR(1,IQ,1) 

IAT=  REGR (1, IQ, 2) 

ENDIF 

315  CONTINUE 

IVAR=0 
I  IN=0 

DO  320  I I -NX, 1,-1 
IVAR=IVAR+ 1 

IF  (MODELS (I AT, 1 1) .EQ. 1)  THEN 
I IN=I IN+ 1 
IBUFF (IIN) =IVAR 
ENDIF 

320  CONTINUE 

SP=RMIN 

WRITE (30 , 1545)  MM.SP, (IBUFF ( I J) ,IJ=1,IIN) 

C 

C  #»*»#»*#*«*#*»#****»#«»**##»»#*#*#*»*»*#»** 

C  »  FIND  THE  VOLUME  REDUCTION  AND  INDICATE  * 

C  *  COVERAGE  » 

C  *#»*«««*»««»#»*#***»»»»##»*#»#»»»**»#»*##»* 

c 

CALL  COVER (VC V , MODELS , KNX , NX , NVAR , I AT , I IN , YBAR , CBAR , 

&  VECMUC , NY , VECMUY , NUMREPS , FF , IH , ICOVER , VOLRED , 

&  VECYBAR, IKNOW,COVCV,VU,DIFF) 

OOHIIIUtflHIHmmHIHimHHtlltl 

»  COVERAGE  AND  VOLUME  REDUCTION  TALLYS  « 

#«#*###»#*#**»*#####*####*##*#»###»###*#*** 

DO  325  IC=  1 ,4 

I CTOT ( I C ) = I CTOT ( I C ) ♦ I  COVER ( IC) 

325  CONTINUE 

DO  330  IC-1,2 

SUMDEV(IC) =SUMDEV (IC) +DIFF 
SUMVU(IC) =SUMVU(IC) +VU 
330  CONTINUE 

ENDIF 

PRINT  « , "THIS  IS  META-EXPERIMENT  •’ ,MM, * ICOVER  ICOVER 
1000  CONTINUE 
C 
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C  *  META  EXPERIMENT  LOOP  ENDS  ON  ABOVE  LINE  « 

C  imumtHimHHimHHmHimHmumi 
C 

DO  1005  12=1 . 2 

SUMDEV(IZ)=SUMDEV(IZ) /FLOAT (META) 

SUMVU ( IZ) =SUMVU ( IZ) /FLOAT (META) 

VR(IZ) =SUMDEV(IZ) /SUMVU(IZ) 

1005  CONTINUE 
C 


1010 

C 


C 


DO  1010  IZ=  1 ,4 

COVERAG ( IZ) = FLOAT ( ICTOT ( IZ) ) /FLOAT (META) 
CONTINUE 


WRITE (30 ,1515) 
WRITE(30, 1550) 
WRITE (30 , 1551) 
WRITE(30, 1555) 
WRITE (30 ,1515) 
WRITE (30 , 1560) 
WRITE(30, 1561) 
WRITE(30 , 1565) 


COVERAG (1) 
VR(  1) 

COVERAG (2) 

COVERAG (3) 
VR  ( 1 ) 

COVERAG (4) 


CLOSE (5) 
CLOSE (10) 
CLOSE (30) 
STOP 


C 

c  MOHIIIIHtllOIOHOHOtOmiltHHO 

C  «  FORMAT  STATEMENTS  (MAIN  PROGRAM)  * 

C  **»»»»»»*»»»*»»»»»»»»»»»»»»*»»«*••»**»»»»*» 

c 

1500  FORMAT (1X.A25, ’META  =  ’,13,’,  NUMREPS  =  ’  ,  13 .  ’ ,  TOTAL  REPS  = 

&I4) 

1505  FORMAT (IX, ’THE  RESPONSE  ARE' , 13X, ’MEAN  ’,14.’  REPS’, 2X, 

&’ STEADY  STATE  MEAN’/) 

1510  FORMAT (2X ,I2,1X,A25,F12.5,4X,F12.5) 

1515  FORMAT ( ’  ’) 

1520  FORMAT ( IX, ’THE  CANDIDATE  CONTROLS  ARE ’, 3X ,’ MEAN  ’,14,’  REPS’, 

&2X, 'STEADY  STATE  MEAN'/) 

1525  FORMAT (/, IX, ’COVARIANCE  MATRIX  OF  CONTROLS  WAS  ESTIMATED’) 

1530  FORMAT ( / , IX , ’ KNOWN  COVARIANCE  MATRIX  OF  CONTROLS  WAS  USED’) 

1535  FORMAT (IX, ’META#’ ,3X, ’  CRITERION  ’, 10X ,’ VARIABLE  SUBSET ’ ) 

1540  FORMAT (1 OX, ’BEST  ’,12,’  REGRESSIONS  WITH  ’,12,’  VARIABLES’//) 

1545  FORMAT (1X,I4,2X,E16.8,10X,30(I2,1X) ) 

1550  FORMAT ( IX, ’CONTRLD  COVERAGE  ON  STEADY  STATE  MEANS  =’,F12.8) 

1551  FORMAT ( IX, ’  AND  VOLUME  REDUCTION  =’,E16.8) 

1555  FORMAT ( IX, ’UNCONTRLD  COVERAGE  ON  STEADY  STATE  MEANS  =’,F12.8) 

1560  FORMAT (IX, ’CONTRLD  COVERAGE  ON  SAMPLE  MEAN  OF  1000  REPS  =’,F12.8) 

1561  FORMAT ( IX, ’  AND  VOLUME  REDUCTION  =’,E16.8) 

1565  FORMAT (IX, ’UNCONTRLD  COVERAGE  ON  SAMPLE  MEAN  OF  1000  REPS= ’ ,F12 . 8) 
1575  FORMAT (IX, ’ *  OF  CANDIDATE  CONTROLS  EXCEED  PROGRAM  LIMIT  OF’, 13) 
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1580  FORMAT (IX, 1  *  OF  RESPONSES  EXCEED  PROGRAM  LIMIT  OF’,13) 

1585  FORMAT (IX, ' »  OF  BEST  REGRESSIONS  TO  KEEP  EXCEED  PROGRAM’) 

1586  FORMAT (IX, ’LIMIT  OF’ ,13) 

1590  FORMAT (IX, ’•  OF  REPLICATIONS  PER  META  EXPERIMENT  EXCEED’) 

1591  FORMAT( IX, ’PROGRAM  LIMIT  OF’,13) 

1595  FORMAT (IX, ’*  OF  META  EXPERIMENTS  EXCEED  PROGRAM  LIMIT  OF’, 14) 
1600  FORMAT (IX, ’PROGRAM  LIMITS  EXCEEDED  BY  ’,11,’  PRIMARY  INPUTS’) 

END 

C 

C  ##*«**#*#****«»*»##*»«#******#**#*#**#*######*#*#***#****#***# 

C  «  END  MAIN  PROGRAM  » 

C  ##»#*«**#*»*##»«**»»#»»*****«##*»#***#***#****#**#»*»*»**##*»* 

c 

C  tmotiooHtimttmiiiHmiMiitHtmiuittiHmoimH 
C  «  SUBROUTINES  * 

C  »«»***»****»#»*»»*»**#***«**#»«*«*«»«**»**»*#*»#*t***«#*»*»*** 

c 

C  *  SUBROUT INF  COVER  * 

C  »  « 

C  »  THIS  SUBROUTINE  DOES  THE  COVERAGE  AND  VOLUME  » 

C  *  REDUCTION  CALCULATIONS  FOR  THE  OPTIMAL  CONTROL  * 

C  »  SUBSET  * 

C  *«««*••««*••«***«***•«********»»*****»***#*»««*»*»*** 

c 

SUBROUT I NE  COVER ( VC V , MODELS , KNX , NX , NVAR , I AT , I I N , YBAR , CBAR 
& , VECMUC , NY , VECMUY , NUMREPS , FF . IH , ICOVER , VOLRED , VECYBAR , IKNOW 
& .COVCV , VU.DIFF) 

C 

INTEGER  Z1 ,  Z2 ,  Z3 ,  Z5,  Z8 

REAL  PI 
C 

PARAMETER  (Zl=8 .  Z2=8,  Z3=Z1*Z2,  Z5=2«*Z1,  Z8= ( (Z3» (Z3+ 1 ) ) /2) ) 

PARAMETER  (PI=3 . 14 15927) 

C 

REAL  VCV (Z8) ,YBAR(Z2) ,CBAR(Z1) .VECMUC(Zl) ,VECMUY(Z2) 

& , FF (0 : Z 1 ) , VECYBAR (Z2) ,VOLRED(2) ,C0VCV(Z1 ,Z1) 

C 

INTEGER  M0DELS(Z5,Z1) ,IH(Z3) ,IC0VER(4) ,IIN 
C 

REAL  SCBAR(Zl) .SVECMU(Zl) ,SUBV(Z8) ,SUBVF (Z3 ,Z3) ,B(Z3) 

& , WKAREA(2«Z3) .BUFFI (Z3 ,Z3) ,BUFF2 (Z2 ,Z1) ,BETA(Z2,Z1) 

&.CDEV1 (1 ,Z1) .CDEV2 (Z1 , 1) ,EXPL(Z2,Z2) ,DEV(Z2,1) , YBHAT (Z2) 
&,BUFF3(Z1 ,Z2) , BUFF4 (Z2 , Z2) ,SYDOTC (Z2 ,Z2) ,HPH(1,1) ,T1 (1 ,Z1) 

& , YMD1 ( 1 ,Z2) , YMD2 (Z2 , 1 ) ,T2(1,Z2) ,0BS (1,1) ,BUFF5 (Z2 ,Z2) 

& , BUFF6 (Z2 ,Z2) ,YMD3(1,Z2) ,YMD4(Z2,1) , OBS2 (1,1) 

& .SYMCOVC ( (Zl« (Zl+1) ) /2) , SUBCOVC ((Zl«(Zl+l))/2) 

& , FULCOVC (Z1 ,Z1) ,GAMM(Z2,Z1) ,EHAT(Z2,Z2) ,BUFF9 (Z2 ,Z2) 

& , CANCORR ( Z2 , Z2 ) ,REIGS(Z2) ,EIGS(2*Z2) .DUMMY (Z2 ,Z2) 

&. WK (Z2) 

C 
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C 


INTEGER  IH2 (ZI) 


COMPLEX  CEIGS (Z2) 

C 

EQUIVALENCE  (EIGS ( 1) .CEIGS ( 1) ) 

C 

C  #########*#**##»#**«##*«########*»#*#»§»*#» 
C  *  INITIALIZE  COVERAGE  AND  VOLUME  » 

C  »  REDUCTION  VECTORS  * 

C 

DO  10  1=1,4 
ICOVER ( I ) =0 
10  CONTINUE 
C 

DO  15  1=1,2 
VOLRED ( I ) =0 . 

15  CONTINUE 
C 

C  #§**##****#*#«»*#»*»*###*«*#»*»#**#»#**#*#* 
C  *  FIND  THE  SUBMATRIX  FOR  THE  SELECTED  • 
C  *  MODEL  * 

C  »***##»#»#»#»»**##*#*#**»*#»*****«*t#»#»*»« 

C 

DO  20  1=1, NVAR 

IF  (I.LE.NX)  THEN 
I H ( I ) =0 
IH2 (I) =0 
ELSE 

I H ( I ) =1 
END  IF 
20  CONTINUE 
C 

IVAR=0 

DO  25  II=NX, 1,-1 
IVAR=IVAR+1 

IF  (MODELS (I AT, II) .EQ. 1)  THEN 
IH ( IVAR) = 1 
IH2 ( IVAR) = 1 
END  IF 
25  CONTINUE 
C 

M1=Z3 

CALL  RLSUBM(VCV,M1 , IH.SUBV.M2) 


»  FIND  THE  SUB VECTOR  (POPULATION  AND  » 

*  SAMPLE)  OF  THE  CONTROL  MEANS  « 


INDEX=0 


72 


o  o  o  o  oooooo  oooooo 


DO  30  11=1, NX 

IF  (IH(II) . EQ. 1)  THEN 
INDEX= INDEX+ 1 
SCBAR ( INDEX) =CBAR (II) 

SVECMU (INDEX) =VECMUC (II) 

END  IF 
30  CONTINUE 
C 

C  a**#*###**#*#**#**#*#*####***##*#*######*#* 

C  »  BUFFER  THE  COVARIANCE  MATRIX  OF  * 

C  *  SELECTED  CONTROLS  AND  RESPONSES  « 

C  in**#****#####***##*#**#**#*****#**#****### 

c 

CALL  VCVTSF (SUBV , M2 , SUBVF , Z3) 

C 

DO  35  1=1 ,M2 
DO  35  J  =  1 , M2 

BUFFI (I , J) = SUBVF (I , J) 

35  CONTINUE 

««***»**«*««««***«*•*•***»««*«»»«•*«««***«* 

*  INVERT  THE  COVARIANCE  SUBMATRIX  OF  » 

»  CONTROLS  « 

»**»»*»*i*#*»»#**»***»***#*#«»«#***»*»»»##* 

CALL  LINV3F (SUBVF ,B, 1 , I IN ,Z3 , D1 , D2 , WKAREA , IER) 

IF  (IER.NE.O)  PRINT  *, ‘I  DIED  BELOW  35  (SUBROUTINE  COVER)’ 

#**#*»#*»**»#*#*«*#*#***»*«#»#*##****»*#*** 

*  BUFFER  THE  CROSS-COVARIANCE  SUBMATRICES  * 

*  OF  SELECTED  CONTROLS  WITH  RESPONSES  » 

HHIIHtXtKHKHIHIHIIIHIXKVHHHIIkDIlK 

DO  40  I=IIN+1 ,M2 
DO  40  J= 1 , I IN 

BUFF2 (I-IIN, J) =BUFF1 (I , J) 

BUFF3 (J , I-IIN) = BUFFI (I , J) 

40  CONTINUE 

#*#**##**»»####* 

«  BUFFER  THE  COVARIANCE  SUBMATRIX  OF  * 

«  RESPONSES  t 

C  HHmmHHHMUmHUHUHMXXIHIM 

C 

DO  45  I=IIN+1 ,M2 
DO  45  J  =  1 1 N ♦ 1 , M2 

BUFF4 (I-IIN, J-IIN) =BUFF 1 (I , J) 

BUFF6 (I-I7N.J-IIN) = BUFFI ( I , J) 

45  CONTINUE 
C 
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C  HHitimitimoommKiHmmimii 
C  *  FIKD  THE  BETA  HAT  MATBIX  (  CONTROL  « 

C  *  COEFFICIENTS  )  OR  THE  GAMMA  HAT  MATRIX  » 

C 

IF  (IKNOW.EQ.O)  THEN 

CALL  VMULFF (BUFF2 , SUBVF , NY , I I N , I IN , Z2 , Z3 , BETA , Z2 , IER) 

ELSE 

CALL  VMULFF (BUFF2 , SUBVF , NY , I IN , I IN , Z2 , Z3 , BETA , Z2 , IER) 

CALL  VCVTFS (COVCV,NX,Zl .SYMCOVC) 

CALL  RLSUBM( SYMCOVC ,NX , IH2 , SUBCOVC , I ORDER) 

CALL  VCVTSF (SUBCOVC , IORDER .FULCOVC ,Z1) 

CALL  LINV3F ( FULCOVC . B , 1 , I IN , Z1 , D 1 , D2 , WKAREA , IER) 

IF  (IER.NE.O)  PRINT  *,*I  DIED  BELOW  45  (SUBROUTINE  COVER)* 
CALL  VMULFF ( BUFF2 , FULCOVC , NY , I IN , I IN , Z2 , Z 1 , GAMM , Z2 , IER) 

END  IF 
C 

C  »  FIND  THE  VECTOR  OF  CORRECTIONS  TO  * 

C  #  CONTROL  Y  BAR  » 

C  #»###«»*#*«#»#»#******»*«***»*#»*»*»*»»»**» 

c 

DO  50  1  =  1, 1 IN 

CDEV1 (1,1) =SCBAR ( I ) -SVECMU ( I ) 

CDEV2 (1,1) =CDEV1 (1,1) 

50  CONTINUE 
C 

IF  (IKNOW.EQ.O)  THEN 

CALL  VMULFF (BETA, CDEV2, NY, I IN, 1 ,Z2 ,Z1 ,DEV,Z2 , IER) 

ELSE 

CALL  VMULFF (GAMM. CDEV2, NY, I IN, 1 ,Z2,Z1 .DEV.Z2.IER) 

ENDIF 

C 

C  ****»***»t«*»«#f**»*t#»»»»***»*»»**»»»»*»«* 

C  *  FIND  THE  CONTROLLED  ESTIMATOR  OF  THE  * 

C  »  MEAN  * 

C  t»*#*««#*t*»*»»**»*«#*#*»*****#«*«»*»*»*»»* 

C 

DO  55  1*1 , NY 

YBHAT ( I ) = YBAR ( I ) -DEV (1,1) 

55  CONTINUE 
C 

C  *  FIND  THE  MATRIX  OF  EXPLAINED  COVARIANCE  » 

C  *  DUE  TO  CONTROL  * 

c  . . nil 

C 

CALL  VMULFF (BETA , BUFF3 , NY , 1 1 N , NY , Z2 , Z 1 , EXPL , Z2 , IER) 


•  FIND  THE  RESIDUAL  COVARIANCE  » 
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c 


C 1 = (FLOAT (NUMREPS- 1 ) /FLOAT (NUMREPS-I IN- I ) ) 


C 

DO  60  I = 1 , NY 
DO  60  J= 1 ,NY 

SYDOTC (I , J) = (BUFF4 (I , J) -EXPL(I ,J) ) *C1 
BUFF5 (I , J)  = SYDOTC (I ,  J) 

60  CONTINUE 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

»  FIND  THE  ESTIMATOR  SIGMA  TILDE  HAT  x 

XXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

IF  (IKNOW.EQ.l)  THEN 

CONST 1= (FLOAT (NUMREPS- 2) ) / (FLOAT (NUMREPS* (NUMREPS- 1 ) ) ) 
CONST2= (FLOAT (IIN+1) ) / (FLOAT (NUMREPS* (NUMREPS- 1 ) ) ) 

DO  65  1=1, NY 
DO  65  J  = 1 , NY 

EHAT(I , J) = (C0NST1»SYD0TC (I , J) ) + (CONST2*BUFF4 (I , J) ) 
BUFF9 ( I , J) =EHAT ( I , J) 

65  CONTINUE 
END  IF 

XXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXX 

X  FIND  THE  INVERSE  RESIDUAL  COVARIANCE  * 
x  MATRIX  * 

XXXXXXXXXXXXXXXXXXXXXKXXXXXKXXXXXXXXXXXXXXX 

IF  (IKNOW.EQ.O)  THEN 

CALL  LINV3F (SYDOTC , B , 1 ,NY,Z2,D1 ,D2 .WKAREA, IER) 

IF  (IER.NE.O)  PRINT  *  ,  ’  I  DIED  BELOW  65  [1]  (SUBR  COVER)’ 
ELSE 

CALL  LINV3F (EHAT ,B , 1 ,NY,Z2,D1 ,D2 .WKAREA, IER) 

IF  (IER.NE.O)  PRINT  »,’I  DIED  BELOW  65  [1 ]  (SUBR  COVER)’ 
END  IF 


x  COMPUTE  THE  DEVIATIONS  FROM  THE 
x  STEADY-STATE  RESPONSE  VECTOR 
*  (BOTH  CASES:  CONTROLLED /UNCONTROLLED) 


DO  70  1=1. NY 

YMD 1(1,1) = YBHAT ( I ) -VECMUY ( I ) 
YMD2 (1,1) =YMD1 (1,1) 

YMD3 (1,1) = YBAR ( I ) -VECMUY ( 1 5 
YMD4 (1,1) =YMD3 (1,1) 

70  CONTINUE 
C 
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c  »»**»«»*»**#*»*»*t»«##***»#**»#*»#*#*»*tHt** 

C  *  COMPUTE  H’H  * 

C  »  (NOTATION  AS  PER  VENXATRAMAN  AND  * 

C  »  WILSON  1986)  » 

C  *#*#»«*»**»»#»#»*»«#»*»»*»»»»»»»#**#****»»* 

c 

IF  (IKNOW.EQ.O)  THEN 

CALL  VMULFF (CDEV1 .SUBVF , 1 , 1 IN , I  IN , 1 ,Z3 ,T1 , 1 . IER) 

CALL  VMULFF (T1.CDEV2 .1,1 IN , 1 , 1 ,Z1 ,HFH. 1 , IER) 

END  IF 
C 

IF  (IKNOW.EQ.O)  THEN 

X= ( 1 . /FLOAT (NUMREPS) ) + ( 1 . /FLOAT (NUMREPS- 1 ) ) *HPH (1,1) 
ELSE 
X=l. 

END  IF 
C 

C  *»»**#»»»»##**»»*#**»#»*»#***#*«*#*«****#*# 

C  »  COMPUTE  THE  RIGHT  HAND  SIDE  FOR  THE  » 

C  »  CONFIDENCE  REGION  AS  PER  RAO  (1967)  « 

C  »»***»*###»*##**»***#»*»*****»**»»**»*#*##* 

c 

C2= (FLOAT ( (NUMREPS- I IN- 1 ) *NY) /FLOAT (NUMREPS- I IN-NY) ) 
F=EXP ( ( 1 . /FLOAT (NY) ) «ALOG(FF (I  IN) ) ) 

RHS=X*C2*F 

C 

C  tt#***«*«f#f#****«»f#****f*##f#«###*«******* 

C  *  COMPUTE  THE  T»*2  STATISTIC  FOR  THE  CASE  » 

C  «  WHERE  CONTROLS  ARE  USED  (STEADY  STATE  * 

C  *  ASSUMED)  # 

C  **««**#«»»##*#»«»*#»*#*»»##*#*******#*»#### 

c 

IF  (IKNOW.EQ.O)  THEN 

CALL  VMULFF (YMD1 , SYDOTC , 1 , NY , NY , 1 , Z2 , T2 , 1 , I ER) 

CALL  VMULFF (T2 , YMD2 , 1 , NY , 1 , 1 , Z2 , OBS , 1 , IER) 

ELSE 

CALL  VMULFF (YMD1 , EHAT , 1 , NY . NY , 1 , Z2 , T2 , 1 , I ER) 

CALL  VMULFF (T2.YMD2.1, NY, 1.1.Z2, OBS, 1, IER) 

ENDIF 


»  INDICATE  COVERAGE  FOR  THE  CASE  WHERE  » 

«  CONTROLS  ARE  USED  (STEADY  STATE  ASSUMED) * 


IF  (OBS (1,1) .LE.RHS)  THEN 
T COVER ( 1) = 1 
ELSE 

ICOVER ( 1 ) =0 
ENDIF 
C 
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C  immmHmoHHHiimHHimmHM 
C  *  COMPUTE  THE  VOLUME  REDUCTION  * 

C 

IF  (IKNOW.EQ.O)  THEN 

CALL  LINV3F (BUFF4 ,B,4,NY1Z21D1 ,D2 .WKAREA, IER) 

IF  (IER.NE.O)  PRINT  «,*I  DIED  BELOW  70  El]  (SUBR  COVER) 
UCDET=D1*2**D2 

CALL  LINV3F (BUFF5 , B , 4 , NY , Z2 , D1 , D2 , WKAREA , IER) 

IF  (IER.NE.O)  PRINT  *,*I  DIED  BELOW  70  [2]  (SUBR  COVER) 
CDET=D1«2**D2 
ELSE 

CALL  LI NV3F ( BUFF4 , B , 4 , NY , Z2 , D 1 , D2 , WKAREA , I ER) 

IF  (IER.NE.O)  PRINT  *,’I  DIED  BELOW  70  13]  (SUBR  COVER) 
UCDET=D1*2**D2 

CALL  LINV3F (BUFF9 , B , 4 , NY ,Z2 , D1 ,D2 .WKAREA, IER) 

IF  (IER.NE.O)  PRINT  *,‘I  DIED  BELOW  70  [4]  (SUBR  COVER) 
CDET=D1«2»*D2 
END  IF 
C 

TERM1 = (CDET/UCDET) »» ( . 5) *X*« (FLOAT (NY) / 2 . ) 

C3=FLOAT ( (NUMREPS-IIN- 1) * (NUMREPS) » (NUMREPS-NY) ) 

C4=FLOAT ( (NUMREPS- 1 IN-NY) « (NUMREPS- 1 ) ) 

TERM2= (C3/C4) »* (FLOAT (NY) /2. ) 

F2=EXP ( ( 1 . /FLOAT (NY) ) *ALOG (FF (0) ) ) 

TERM3= (F/F2) »» (FLOAT (NY) /2 . ) 

VOLRED ( 1 )  =  ( 1 . - (TERM1»TERM2«TERM3) ) *  100 . 

C 

C  HMHHMOilltOtiUimmtlMlIHHOt* 

C  »  COMPUTE  THE  ACTUAL  VOLUME  OF  THE  » 

C  »  CONTROLLED  ELLIPSOID  t 

C 

POVER2=FLOAT(NY)/2. 

CC 1=1. /FLOAT (NY) 

CC2= (2 . «PI*«P0VER2) /GAMMA (P0VER2) 

CC3=FL0AT(NY« (NUMREPS- I IN- 1) ) 

CC4=FL0AT (NUMREPS- I IN-NY) 

CC3= (CC3/CC4) *«P0VER2 
CC4=SQRT (FF (I IN) ) 

VC=CC 1 *CC2*CC3*CC4«SQRT (CDET) « (X»»POVER2) 


*  COMPUTE  THE  ACTUAL  VOLUME  OF  THE  » 

»  UNCONTROLLED  ELLIPSOID  « 


CC3=FL0AT(NY» (NUMREPS- 1 ) ) 

CC4=FLOAT (NUMREPS* (NUMREPS-NY) ) 

CC3=CC3/CC4 

VU=CC1»CC2*SQRT (UCDET) » (CC3«FF ( I IN) ) **POVER2 
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c  **»»#»»*#«»*»*»**»##*»«*##«***»»**#»#»»»#»» 

C  »  COMPUTE  THE  DIFFERENCE  (DIFF)  « 

C  #**«*»»#»»»*»«*#»#»»****»***»»»*»»*##»*»#** 

c 

DIFF=VU*VC 

C 

C  KiHiHOHOiotiiMiHimtmmmMiti 
C  »  COMPUTE  THE  T»«2  STATISTIC  FOR  THE  CASE  » 

C  *  WHERE  NO  CONTROLS  ARE  USED  « 

C  #*«##*»**»*#*******»**###*##**»#*»*»#»*#«** 

c 

CALL  LINV3F (BUFF6 , B , 1 , NY ,Z2 , D1 ,D2 , WKAREA, IER) 

IF  (IER.NE.O)  PRINT  »,'I  DIED  BELOW  70  [5 ]  (SUBR  COVER)' 
CALL  VMULFF (YMD3 , BUFF6 , 1 . NY , NY , I , Z2 , T2 , 1 , IER) 

CALL  VMULFF (T2 ,YMD4 , 1 , NY ,1,1 ,Z2 ,0BS2 . 1 , IER) 

C 

C  «#»«»*#****#««**«f******#***#*»*##*****##»* 

C  »  COMPUTE  THE  RIGHT  HAND  SIDE  FOR  THE  « 

C  »  CONFIDENCE  REGION  * 

C  *»»»*##**»#»»#***»»******#**»*»*»»*#*»*#*** 

c 

C5= (FLOAT ( (NUMREPS-1) »NY) /FLOAT ( (NUMREPS- NY) »NUMREPS) ) 
RHS2=EXP ( ( 1 . /FLOAT (NY) ) #ALOG (FF (0) ) ) *C5 
C 

C  *«***»***«««»«»«**##*»»«#*«»»«»«**»#*•*»**» 

C  *  INDICATE  COVERAGE  FOR  THE  CASE  WHERE  * 

C  *  NO  CONTROLS  ARE  USED  (STEADY  STATE  » 

C  »  ASSUMED)  « 

C  »*****#****»**»#»##»#**»*tHt»**#****«****t#« 

c 

IF  (OBS2 ( 1 ,1) .LE.RHS2)  THEN 
ICOVER(2) = 1 
ELSE 

ICOVER(2) =0 
END  IF 
C 

C  *»*»»»»*»»»**#*»*»»»»#*«*#»*«***»#*#»**#### 

C  *  THE  REMAINING  ANALYSIS  DUPLICATES  THE  « 

C  *  ABOVE,  EXCEPT  THAT  THE  GRAND  MEAN  OF  * 

C  »  1000  RESPONSES  IS  USED  » 

C 

c  HimmHmmimmHHOHOHiHHH 
C  »  RECOMPUTE  DEVIATIONS  » 

C 

DO  75  1*1, NY 

YMD1 (1,1) *YBHAT(I) -VECYBAR(I) 

YMD2 (1,1) * YMD 1(1,1) 

YMD3 (1,1) =YBAR ( I ) - VECYBAR ( I ) 

YMD4 (1,1) =YMD3 (1,1) 

75  CONTINUE 
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*  COMPUTE  THE  T»«2  STATISTIC  FOR  THE  » 

«  CASE  WHERE  CONTROLS  ARE  USED  » 

*  (GRAND  MEAN  USED)  » 

oioimmmiHHmmmiHmimiiH 

IF  (IKNOW.EQ.O)  THEN 

CALL  VMULFF (YMD1 , SYDOTC , 1 , NY , NY , 1 , Z2 , T2 , 1 , I ER) 
CALL  VMULFF (T2 , YMD2 , 1 ,NY, I , 1 ,Z2 ,OBS , 1 , IER) 

ELSE 

CALL  VMULFF (YMD1 , EHAT , 1 , NY , NY , 1 , Z2 , T2 , 1 , IER) 
CALL  VMULFF (T2 , YMD2 , I , NY , 1 , 1 ,Z2 , OBS , 1 , IER) 

END  IF 

»  INDICATE  COVERAGE  FOR  THE  CASE  WHERE  » 

«  CONTROLS  ARE  USED  « 

«***•»*«»**«*«•«****«*«««*••«•**«•**«*»«««« 

IF  (OBS(l . 1) .LE.RHS)  THEN 
ICOVER (3) = 1 
ELSE 

ICOVER (3) =0 
END  IF 

•«••«•**«*««•**«••*•*«•*«•••««•««««•«*««•«* 

»  COMPUTE  THE  T*»2  STATISTIC,  FOR  THE  » 

*  CASE  WHERE  NO  CONTROLS  ARE  USED  » 

*  (GRAND  MEAN  USED)  * 

*ft**»*t*«***»«*«*it«***K*ft**ftftft«*««»******«* 

CALL  VMULFF (YMD3 , BUFF6 , 1 , NY , NY , 1 , Z2 , T2 , 1 , I ER) 

CALL  VMULFF (T2 , YMD4 , 1 , NY , 1 , 1 , Z2 , 0BS2 , 1 , IER) 

»  INDICATE  COVERAGE,  FOR  THE  CASE  WHERE  * 

*  NO  CONTROLS  ARE  USED  » 

IF  C0BS2 (1,1) .LE.RHS2)  THEN 
ICOVER (4) = 1 
ELSE 

ICOVER (4) =0 
END  IF 


*  THIS  SECTION  COMPUTES  THE  CANONICAL 

*  CORRELATIONS  FOR  THE  SUBSET  MODELS  AND 
«  THE  FEASIBILTY  BOUND  FOR  USING  THE 

«  KNOWN  COVARIANCE  MATRIX  OF  CONTROLS 
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c 

IF  ( I KNOW. EQ . 1 )  THEN 

CALL  VMULFF ( BUFF6 , EXPL , NY , NY , NY , Z2 , Z2 , CANCORR , Z2 , I ER) 

CALL  El GRF (CANCORR , NY . Z2 , 0 , E I GS , DUMMY , Z2 , WK , I ER ) 

C 

IC0UNT=0 
DO  80  I = 1 , NY 
DO  80  J= 1 , 2 

I COUNT = I COUNT +1 

IF  (J .EQ. 1)  REIGS ( I ) =SQRT (EIGS { ICOUNT) ) 

80  CONTINUE 
C 

CTOP=FLOAT( (NUMREPS+IIN-1) » (NUMREPS-IIN-2) ) / 

&  FLOAT ( (NUMREPS-1) * (NUMREPS-2) ) 

CBOT= CTOP  * ( FLOAT ( NUMREPS - 2 ) / FLOAT ( NUMREPS ♦ 1 1 N- 1 ) ) 

BOUND=SQRT ( (CTOP- 1 . ) / (CBOT-  1 . ) ) 

PRINT  *, ‘CANONICAL  CORRELATIONS  * .REIGS, *  BOUND  ‘ .BOUND 
PRINT  « .EIGS 
ENDIF 
C 

RETURN 

END 

C 

C  *#«*#***#*»«**»*#**#*#*****#»*»****»*»#»»*»#**##**### 

C  *  SUBROUTINE  COVKNOW  * 

C  *  * 

C  *  THIS  SUBROUTINE  RETURNS  THE  GENERALIZED  VARIANCE  » 

C  *  OF  SIGMA  TILDE  HAT  * 

C  »*##*»»»»***##*#t»#»»«»**»»*#*#*»*#*#*##**###*#«***** 

c 

SUBROUT I NE  COVKNOW ( RSS , NY . Z2 , FULL , NVAR , Z3 , TARGET , DUM . NUMREPS , MV 
&  ,DET) 

C 

INTEGER  Z2.Z3, NY, NVAR, NUMREPS 

REAL  RSS (Z2 ,Z2) ,FULL(Z3,Z3) , TARGET (Z2 ,Z2) ,DUM(Z2) 

C 

Cl  = (FLOAT (NUMREPS-2) /FLOAT (NUMREPS* ( NUMREPS -MV- 1 ) ) ) 

C2= (FLOAT (MV+ 1 ) /FLOAT (NUMREPS* (NUMREPS-1 ) ) ) 

NX=NVAR-NY 

C 

DO  10  1=1, NY 
DO  10  J  = 1 , NY 

TARGET (I , J) = (Cl »RSS (I , J) ) + (C2*FULL (NX+I ,NX*J) ) 

10  CONTINUE 
C 

CALL  LINV3F (TARGET , DUM, 4 ,NY,Z2 ,D1 , D2 , WKAREA , IER) 

IF  (IER.NE.O)  PRINT  *,‘I  DIED  BELOW  10  (SUBR  COVKNOW)' 
DET=D1»2**D2 
C 

RETURN 

END 

C 
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C  HomimiotmoHHtHtmmiomiiMMttMHi 


C  *  SUBROUTINE  FTABL  » 
C  »  « 
C  »  THIS  SUBROUTINE  COMPUTES  AN  F  TABLE  * 
C  »  (TO  THE  POWER  P)  « 


C  '••••*****»***««t«*****««««*«***«**«**»*»«««t«*««««ft* 

c 

SUBROUTINE  FTABL(FF,NX,Z1) 

C 

COMMON  /BLK1/  SIG.KK, IQQ, IP 
C 

INTEGER  Z1 ,KK,IQQ,IP,NSIG,NROOT,ITMAX 
C 

REAL  ROOT ( 1 ) ,LAST,FF(0:Z1) ,EPS,FP 
C 

EXTERNAL  F 
C 

EPS= . 001 
NSIG=5 
NROOT= 1 
ITMAX= 1000 
LAST=3 . 

C 

DO  10  IQQ=0 ,NX 
ROOT ( 1 ) =LAST 

15  CALL  ZREAL2 (F ,EPS , EPS ,EPS ,NSIG , NROOT .ROOT , ITMAX, IER) 

IF  (IER.EQ.33)  THEN 
ROOT ( 1 ) =LAST+ 1 . 

IER=0 

WRITE (6 , 1535) 

WRITE (6 , *)  •' 

GOTO  15 
END  IF 

LAST=ROOT ( 1) 

FP=R00T(1)*#IP 
FF (IQQ) =FP 
10  CONTINUE 
C 

RETURN 

1535  FORMAT (IX, ’IGNORE  LAST  IER=33  WARNING  ---  REINITIALIZING') 
END 
C 

C  »  SUBROUTINE  GAUSS  » 

C  *  t 

C  *  THIS  SUBROUTINE  PERFORMS  THE  PIVOTS  FOR  VARIABLE  * 

C  *  INTRODUCTION  INTO  REGRESSION  MODELS:  * 

C  »  FURNIVAL  AND  WILSON  1974  » 

C  . . . 

C 

SUBROUTINE  GAUSS ( IB , IS , IP , A ,KP .MAX ,Z3) 

C 
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INTEGER  I B , I S , I P , KP , MAX , Z3 
C 

REAL  A (MAX ,Z3 ,Z3) 

iioiifHiiHHmmxixifmtoioHito 

«  TOLERANCE  CHECK  ON  PIVOTS  # 

x«xxx«xxx«x««xx***«txxxx«xx 

LB=IP+ 1 

IF  (A ( IB . IP , IP) .LT. .01)  THEN 
DO  10  L=LB,KP 

A(ISIIP,L)=A(IB,IP,L) 

DO  10  M=L ,KP 

A(ISIL,M)=A(IB,L,M) 

10  CONTINUE 

ELSE 

DO  15  L=LB ,KP 

A(IS,IP,L)=A(IB1IP,L)/A(IB,IP,IP) 

DO  15  M=L,KP 

A(IS,L,M)=A(IB,L,M)-A(IB,IP,M)»A(IS,IP.L) 
15  CONTINUE 

END  IF 
C 

RETURN 

END 

C 

C  xxtxtxxxtxxittxixxxxtxixtxxxtxxx 


C  *  SUBROUTINE  KEEP IT  * 
C  *  * 
C  *  THIS  SUBROUTINE  FINDS  THE  MODEL  OF  A  CANDIDATE  » 
C  •  REGRESSION  « 


C  ItXIXItXtXXIXXXIXiXXXXXtXtXXXtIXtl 

c 

SUBROUTINE  KEEPIT (NUMREG ,NK , NX .MODELS ,Z1 ,Z3,Z5) 

C 

INTEGER  Z1 ,Z3,Z5, NX .NUMREG 
INTEGER  NK(Z3) , MODELS (Z5 ,Z1 ) 

C 

DO  10  1=1, NX 

MODELS (NUMREG , I) =NK ( I ) 

10  CONTINUE 
C 

RETURN 

END 

C 

c  XX«XIXXXXXXXXXXIXXX>XXXXXX«XX 

C  »  THE  FOLLOWING  SUBROUTINES,  USED  IN  THIS  « 

C  »  PROGRAM,  ARE  IMSL  ROUTINES:  « 

C  »  # 

C  »  -  BECOVM  « 

C  *  COMPUTES  MEANS  AND  VARIANCE-COVARIANCE  » 

C  »  MATRIX  « 
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C  *  -  EIGBF  « 

C  »  COMPUTES  EIGENVALUES  AND  (OPTIONALLY)  » 

C  *  EIGENVECTORS  FOR  A  REAL  GENERAL  MATRIX  * 

C  *  IN  FULL  STORAGE  MODE  « 

C  *  -  LINV3F  » 

C  »  COMPUTES  IN-PLACE  INVERSE,  EQUATION  * 

C  *  SOLUTION,  AND/OR  DETERMINANT  EVALUATION  * 

C  *  IN  FULL  STORAGE  MODE  » 

C  *  -  MDFD  * 

C  «  COMPUTES  F  PROBABILITY  DISTRIBUTION  * 

C  *  FUNCTION  » 

C  *  -  RLSUBM  * 

C  *  PERFORMS  RETRIEVAL  OF  A  SYMMETRIC  * 

C  «  SUBMATRIX  FROM  A  MATRIX  STORED  IN  « 

C  «  SYMMETRIC  MODE  * 

C  *  -  VCVTFS  » 

C  *  PERFORMS  STORAGE  MODE  CONVERSION  OF  « 

C  »  MATRICES  (FULL  TO  SYMMETRIC)  * 

C  *  -  VCVTSF  * 

C  *  PERFORMS  STORAGE  MODE  CONVERSION  OF  * 

C  »  MATRICES  (SYMMETRIC  TO  FULL)  * 

C  »  -  VMULFF  * 

C  «  PERFORMS  MATRIX  MULTIPLICATION  (FULL  » 

C  «  STORAGE  MODE)  * 

C  *  -  ZREAL2  » 

C  *  COMPUTES  THE  REAL  ZEROS  OF  A  REAL  « 

C  *  FUNCTION  -  TO  BE  USED  WHEN  INITIAL  » 

C  »  GUESSES  ARE  GOOD  « 

C  *»»*»**##*»«#»»***#*»**#**•#*##*»»*»»****»#**»»##»»** 

c 

C  »*»*****##»####»»***#******»**»*»#**#***»#«***t**#t*»*«#****»* 


c 

c 

c 

c 

c 

c 

c 


*  THE  FOLLOWING  FUNCTIONS  ARE  USED  TO  COMPUTE  THE 

*  SELECTION  CRITERION 


«*«»«•*«**«**»#••«««•§**«*««»*««•*«*«****•* 
»  FUNCTION  Cl  « 


REAL  FUNCTION  Cl (K, IQ, IP) 

PROD= 1 . 

DO  10  1  =  1, IP 
ITOP= (K-IQ-I) 

IBOT= (K-IQ- 1 ) *K 
TERM= FLOAT ( ITOP) /FLOAT ( IBOT) 
PROD=PROD*TERM 
10  CONTINUE 
C1=PR0D 
RETURN 
END 


C 


» 

* 
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c  #»*»»»»****«***#****#»»***#*#*##**#*»****»* 

C  *  FUNCTION  C2  * 

C  #**»##*»*»*##«#***#*******»#«*#»*»#*###»*** 

c 

REAL  FUNCTION  C2(K,IQ,IP) 

SUM=0. 

P 1  =  1 . 

P2=  1 . 

DO  10  J=0,IP 

ILEFT=JCOMB (IP , J) 

IF  (J.NE.O)  THEN 

P1=P1* (IQ+2» (J-l) ) 
P2=P2*(K-IC-(2»J)) 

RNEXT=P1/P2 
ELSE 

RNEXT= 1 . 

END  IF 

TERM= FLOAT (I LEFT) *RNEXT 
SUM=SUM+TERM 
10  CONTINUE 
C2=SUM 
RETURN 
END 

**•*«•#«•**«««•••*««*•«•«««•«*«««•**•*•«»•• 

*  FUNCTION  C3  i 

»»»«•«**•**»«*»«*««»*«»*«««*«*««««»«*««*«** 

REAL  FUNCTION  C3(K,IQ,IP) 

C3=C1 (K,IQtIP)*C2(K,IQ,IP) 

RETURN 
END 

«*«**•«**«•»«*«*•**««•««**»«««•»•«**»*•*«** 

*  FUNCTION  C4  * 

»t#»**#»»f#***#»#*****I#**#*K##»*#tft*###»# 

REAL  FUNCTION  C4 (K . IQ . IP) 

PROD= 1 . 

DO  10  1  =  1, IP 

TOP=FLOAT (K- IQ- I ) 

BOT=FLOAT(K-IQ- 1 ) 

PR0D=PR0D* (TOP/BOT) 

10  CONTINUE 
C4=PROD 
RETURN 
END 


»  FUNCTION  C5 
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REAL  FUNCTION  C5(KK,IQ,IP) 

PR0D= 1 . 

DO  10  I=1,IP 

PR0D=PR0D* (FLOAT (KK-IQ- 1 ) /FLOAT (KK- IQ-I ) ) 
10  CONTINUE 
C5=PR0D 
RETURN 
END 

»**•«*•««*««««*«««**««««»«««•*»»•••«•«*••«« 

*  FUNCTION  CFRONT  » 

•♦•a##***#*#*#*#*#*###*#***#****##*##**##*# 

REAL  FUNCTION  CFRONT (K , IQ , IP) 

TOP5FLOAT (K-IQ- 1 ) 

EOT=FLOAT (K- IQ- IP) 

CFRONT5 (TOP/BOT) **IP 

RETURN 

END 

##*#*«**#***#*#«****#**»»#»*»*»**##»»***#*# 

*  FUNCTION  F  « 

*»*«««**tt*t***«*»*«**»«ft«»«tt**««««*****»*» 


REAL  FUNCTION  F(Z) 

COMMON  /BLK1/  SIG, KK , IQQ , IP 
N1  =  IP 

N2=KK-IQQ-IP 

CALL  MDFD (Z , N1 ,N2,P,IER) 

F=SIG-P 

RETURN 

END 


«  FUNCTION  JCOMB  » 


INTEGER  FUNCTION  JCOMB (N,M) 
ITOP=NFACT (N) 

IB0T=NFACT (N-M) * NFACT (M) 

JC0MB=IT0P/IB0T 

RETURN 

END 


*  FUNCTION  NFACT 


* 


INTEGER  FUNCTION  NFACT (M) 
IF  (M.EQ.O)  THEN 
NFACT5 1 
RETURN 
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END  IF 
IP=M 

ILOOP=M- I 

DO  10  I  =  ILOOP  ,2,-1 
IP=IP»I 
10  CONTINUE 
NFACT=IP 
RETURN 
END 
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I.  INTRODUCTION 


Background  And  Purpose 

When  dealing  with  computer  simulations  it  is  typically  desirable  to 
have  a  general  understanding  of  how  the  simulation  inputs  will  affect 
the  final  results.  It  is  also  desirable  to  be  able  to  accurately 
estimate  the  expected  simulation  response.  Furthermore,  if  the 
estimation  of  the  response  can  be  achieved  with  a  subset  of  the 
simulation  inputs  (variables) ,  a  variance  reduction  on  the  estimator  of 
the  mean  can  also  be  realized.  One  way  of  achieving  these  goals  is 
through  the  identification  of  a  good  subset  of  control  variates. 

Control  variates,  also  known  as  control  variables,  are  variables  which 
have  a  significant  covariance  with  the  response(s)  of  interest. 

The  development  of  a  quick  and  easy  method  for  identifying  the 
subset  of  significant  control  variates  in  a  simulation  model  would 
greatly  decrease  the  time  and  effort  required  to  gain  insights  into  the 
simulation.  Identifying  the  significant  control  variates  for  a 
simulation  model  can  also  enhance  the  process  of  preparing  and 
implementing  an  experimental  design.  It  would  eliminate  the  guesswork 
in  determining  which  variables  to  concentrate  on  in  a  subsequent 
experimental  design.  This  could  also  save  computer  time  by  identifying 
a  subset  of  the  available  control  variates  to  work  with,  since  the 

Jt 

standard  experimental  design  requires  2  simulation  runs  to  acquire 
data,  where  k  is  the  number  of  variables  being  tested. 

The  corresponding  purpose  of  the  Variable  Subset  Selection  Program 
(VSSP)  is  to  provide  a  means  to  identify  the  significant  control 
variables,  of  a  simulation,  by  evaluating  the  simulation  output. 


Program  Description 

The  Variable  Subset  Selection  Program  (VSSP)  identifies  the 
significant  control  variables  using  the  ’Best  Controls'  (i.e.  BCp) 
criterion  developed  by  Bauer  and  Wilson  (1990).  Initially,  the  best  or 
near-best  subset  of  variables,  depending  on  the  evaluation  procedure 
desired,  is  selected  for  each  subset  size  from  1  to  MX  (the  total  number 
of  control  variables  in  the  full  set) .  The  initial  selection  among 
subsets  of  the  same  size  is  based  on  the  RSS  (Residual  Sums  of  Squares) 
of  each  subset. 

When  these  subsets  are  selected,  the  corresponding  BCp  criterion 
value  is  calculated  and  the  subset  with  the  best  criterion  value  is 
selected.  The  major  advantage  of  using  the  BCp  criterion  is  that  it 
takes  the  number  of  variables  in  the  subset  into  account  in  determining 
the  criterion  value.  Unlike  other  more  common  selection  criterions,  the 
BCp  criterion  does  not  automatically  select  the  subset  with  the  most 
variables . 

Upon  selection  of  the  'best'  variable  subset,  the  corresponding 
coverage  and  volume  reduction  of  the  confidence  region  is  determined. 

The  values  are  summed  up  over  each  meta-experiment  and  averaged  at  the 
conclusion  of  the  program  rum.  These  values  along  with  the  "best" 
criterion  value  and  variable  subset  for  each  meta-experiment  are  written 
to  an  output  file. 

As  mentioned  earlier,  the  resulting  variable  subset  may  be  the  best 
or  near-best  one  possible,  depending  on  the  selection  procedure  desired. 
The  user  may  choose  from  two  selection  procedures:  1)  Enumerated 
Subsets,  or  2)  Stepwise  (Forward  Selection).  The  enumerated  subsets 
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procedure  evaluates  all  possible  combinations  of  variables,  for  each 
subset  size.  This  ensures  that  the  best  subset  will  be  selected,  based 


on  the  BCp  criterion. 

The  stepwise  (forward  selection)  procedure  does  not  evaluate  all 
possible  subsets,  except  for  the  one-variable  subset  case.  Initially, 
it  evaluates  all  one-v'riable  subsets  then  selects  the  best  one.  After 
selecting  the  best  single  variable  the  procedure  then  evaluates  only 
those  two-variable  subsets  containing  that  one  variable;  selecting  the 
best  of  these  two-variable  subsets.  Then  the  procedure  only  evaluates 

the  three-variable  subsets  containing  the  variables  of  the  best  two- 
variable  subset;  selecting  the  best  out  of  those.  This  process 
continues,  building  on  the  last  subset  selected,  until  all  variables  are 
in  the  model.  The  disadvantage  of  this  procedure  is  it  ignores  all 
subsets  which  do  not  include  the  previously  selected  variables;  so 
better  variable  subsets  may  be  missed.  However,  the  advantage  of  this 
procedure  is  that  an  efficient  implementation  will  select  the  near-best 
variable  subset  in  much  less  time  than  the  enumerated  subsets  procedure. 

In  addition,  the  VSSP  provides  the  user  with  several  more  options. 
The  user  may  specify  the  following: 

-  Whether  the  covariance  matrix  of  control  variables  should  be 
estimated,  based  on  the  data,  or  provide  the  covariance  matrix, 
if  known. 

-  Whether  the  program  should  provide  the  single  best  variable 

subset  or  the  'M‘  best  X-variable  subsets  (X  =  1,2 . NX)  for 

each  meta-experiment, 

-  Whether  or  not  program  input  will  be  provided  by  datafile  or 
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input  manually. 


Brief  Overview  Of  Contents 

Following  sections  of  this  manual  are  designed  to  help  the  user 
understand  the  operation  and  constructs  used  in  the  VSSP.  Section  II 
provides  information  on  how  to  run  the  VSSP  and  the  options  available. 
Section  III  lists  all  the  variables  used  in  the  program  and  explains 
the  purpose  of  each.  Section  IV  provides  the  parameters  associated 
with  each  program  variable  (i.e.  variable  type,  precision  (single  or 
double) ,  whether  or  not  it  is  an  array  and  corresponding  size)  and  any 
applicable  comments.  Section  V  provides  a  listing  of  all  the 
subroutines  used  in  the  program  and  a  brief  description  of  each.  And, 
Section  VI  provides  a  listing  of  all  the  functions  used  in  the  program 
and  a  brief  description  of  each.  In  addition,  several  appendixes  are 
provided  to  supplement  the  information  contained  in  the  various  manual 
sections . 
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II.  HOW  TO  USE  THE  PROGRAM  AND  INTERPRET  THE  OUTPUT 

Program  Operation  And  Required  Data 

When  running  the  Variable  Subset  Selection  Program  (VSSP)  the  user 
will  be  prompted  for  the  following  data  and  information.  The  minimum 
required  data  is  a  file  containing  simulation  output  corresponding  to 
the  control  variables  and  responses  to  be  evaluated  by  the  program.  The 
data  must  be  arranged  with  the  variable  output  values  before  the 
response  values.  Other  data  which  the  user  may  input  to  the  program,  if 
known,  is  the  covariance  matrix  between  the  control  variables. 

The  other  information,  required  by  the  program  to  operate,  is 
identified  in  the  following  list: 

-  Whether  program  data  and  information  will  be  input  by  datafile  or 
manually.  If  the  datafile  option  is  chosen,  the  user  will  be 
prompted  for  the  datafile  name  so  program  data  and  information 
can  be  read  in  by  the  program.  Regardless,  the  following 
information  is  still  required. 

-  Number  of  control  variables. 

-  Number  of  responses. 

-  Number  of  best  regressions  to  keep.  Input  a  one  (1)  if  only  the 
best  variable  subset  is  desired. 

-  Number  of  data  replications  per  meta-experiment. 

-  Number  of  meta-experiments. 

-  Whether  or  not  the  covariance  matrix  of  controls  is  known  or  is 
to  be  estimated  from  the  the  data. 

-  Which  evaluation  procedure  you  desire,  enumerated  subsets  or 
stepwise. 
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-  Level  of  significance  to  use  in  deriving  F-values. 

-  Known  or  estimated  means  of  the  control  variables. 

-  Sample  means  of  the  control  variables. 

-  Estimated  means  of  the  responses. 

-  Sample  means  of  the  responses. 

-  A  title  to  write  to  the  program  summary  file.  For  this  and  any 
other  input  where  character  data  is  necessary,  use  an  underscore 

instead  of  spaces.  Also,  the  input  must  begin  with  a 
character,  not  a  number. 

-  A  name  for  each  control  variable. 

-  A  name  for  each  response. 

-  A  name  for  the  program  summary  datafile. 

If  the  covariance  matrix  of  controls  is  known  and  is  to  be 
provided,  the  program  will  ask  if  this  information  will  be  entered  by 
datafile  or  manually.  If  the  datafile  option  is  selected,  then  the  user 
will  be  prompted  to  provide  the  name  of  the  datafile  where  this 
information  is  contained.  When  the  covariance  matrix  is  to  be  provided 
and  the  program  input  is  by  datafile  then  the  covariance  data  may  be 
provided  by  a  separate  datafile  or  contained  in  the  initial  datafile 
(See  Appendix  B  for  further  detail  on  datafile  format). 

Finally,  there  are  two  items  to  remember  in  executing  the  program 
or  it  may  not  work  properly.  First,  the  number  of  replications  per 
meta-experiment  must  be  greater  than  the  number  of  control  variables 
plus  two  (i.e.  NX+2) .  And  second,  the  datafile  containing  the  control 
and  response  data  must  contain,  as  a  minimum,  a  number  of  data  sets 
equal  to  the  total  replications  (i.e.  'Number  of  meta-experiments’  times 
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’Number  of  replications  per  meta-experiment) . 

Interpreting  The  Output 

After  the  program  has  completed  evaluation  of  the  data,  all  results 
are  output  to  a  summary  file.  The  summary  file  will  provide  the 
following  information  (starting  at  top-left  of  file  and  moving  to  the 
right) : 

-  Evaluation  title,  input  by  user  during  program  data  input  phase. 

-  Number  of  meta-experiments  performed. 

-  Number  of  data  replications  per  meta-experiment. 

-  Total  number  of  data  replications  evaluated  by  the  program.  This 
is  equal  tc  ’Number  of  meta-experiments’  times  ’Number  of 
replications  per  meta-experiment’. 

-  Summary  of  response  data  including  response  number,  name 
designation  of  response,  response  mean  over  total  number  of 
replications,  and  estimated  steady  state  mean.  This  information 
is  provided  by  user. 

-  Summary  of  control  variable  data.  Covers  same  information  as  for 
responses.  This  information  is  provided  by  the  user. 

-  A  statement  whether  the  covariance  matrix  of  controls  was 
estimated  or  known. 

-  The  meta-experiment  number,  best  BCp  criterion  value  found  for 
the  meta-experiment,  and  the  corresponding  ‘best’  variable 
subset. 

-  Two  sets  of  coverage  and  volume  reduction  data  averaged  over  all 
the  meta-experiments.  The  first  set  is  based  on  the  steady  state 
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means  and  the  second  set  is  based  on  the  sample  means. 

The  following  should  be  noted  in  regards  to  the  coverage  and  volume 
reduction  summary  data.  First,  if  the  steady  state  and  sample  means 
input  to  the  program  are  the  same,  there  will  be  no  difference  in  the 
values  of  either  set.  Also,  each  set  of  coverage  and  volume  reduction 
will  contain  two  values  for  coverage.  The  primary  value  of  interest  is 
the  controlled  coverage.  This  corresponds  to  data  coverage  achieved  by 
using  the  selected  variable  subsets. 
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III.  VARIABLE  DICTIONARY 


A  =  Storage  for  full  covariance  matrix  and  subsequent  pivots. 

ANSWER  -  Character  variable  used  to  input  answers  to  initial  data 
input  questions. 

B  =  Dummy  array  for  IMSL  Subroutine  LINV3F.  _. 

BETA  =  Equivalent  to  BUFF2  »  SUBV  =  &  =  S„C«SCC 

BIG  =  Equivalent  to  [ (NUMREPS- 1 ) / (NUMREPS -NX- 2) )**NY;  partial  value 

used  in  computing  TWO. 

BOT  =  Denominater  portion  of  a  value  used  to  compute  C4  and  CFRONT, 

in  their  respective  functions. 

BOUND  =  Feasibility  bound  when  using  known  covariance  matrix  of 
controls.  Used  in  Subroutine  COVER. 

BUFF  =  Buffer  used  in  book  keeping  of  REGR. 

BUFFI  =  Buffer  for  S  (Full  Covariance  Matrix),  [See  Note  1]. 

BUFF2  =  Buffer  for  Syc,  (See  Note  1]. 

BUFF3  =  Buffer  for  sL .  (See  N"te  1]. 

BUFF4  =  Buffer  for  Syy ,  (See  Note  1]. 

BUFF5  =  Buffer  of  SYDOTC. 

BUFF6  =  Buffer  of  Syy ,  (See  Note  1). 

BUFF9  =  Buffer  of  Sigma  Tilde  Hat. 

Cl  =  Equivalent  to  (NUMREPS- 2) / (NUMREPS- I IN- 1 ) ;  partial  value  of 

residual  covariance.  Used  in  Subroutine  COVER.  Also  used  in 
Subroutine  COVKNOW  to  find  generalized  variance  of  Sigma 
Tilde  Hat. 

C2  *  Equivalent  tc  (MV* 1 )/ (NUMREPS (NUMREPS- 1 )) ;  partial  value  of 

generalized  variance  of  Sigma  Tilde  Hat.  Used  in  Subroutine 
COVKNOW.  Also  used  in  Subroutine  COVER  in  computing  the 
right  hand  side  of  confidence  region  when  controls  are  used. 

C3  =  Equivalent  to  (NUMREPS- I IN- 1 )» (NUMREPS-NY) ;  partial  value 

used  in  computing  volume  reduction  in  Subroutine  COVER. 

C4  =  Equivalent  to  (NUMREPS-I IN-NY) * (NUMREPS- 1) ;  partial  value 

used  in  computing  volume  reduction  in  Subroutine  COVER. 

C5  =  Equivalent  to  (NY (NUMREPS- 1 ))/ (NUMREPS (NUMREPS-NY) ) ;  partial 

value  used  in  computing  right  hand  side  for  the  confidence 
region  where  no  controls  are  used.  Used  in  Subroutine  COVER. 

CANCORR  =  Canonical  Correlations  for  variable  subset  models  when  using 
a  known  covariance  matrix  of  controls.  Used  in  Subroutine 
COVER. 

CBAR  =  Sample  mean  vector  for  controls. 

CBOT  =  Denominator  of  equation  used  in  calculating  bound  when  the 
covariance  matrix  of  controls  is  known.  Used  in  MAIN 
program. 

CC1  *  Partial  value,  used  to  compute  actual  volume  of  the 

controlled  and  uncontrolled  ellipsoids.  Used  in  Subroutine 
COVER. 

CC2  =  Partial  value,  used  to  compute  actual  volume  of  the 

controlled  and  uncontrolled  ellipsoids.  Used  in  Subroutine 
COVER. 


CC3  -  Partial  value,  used  to  compute  actual  volume  of  the 

controlled  and  uncontrolled  ellipsoids.  Used  in  Subroutine 
COVER. 

CC4  =  Partial  value,  used  to  compute  actual  volume  of  the 

controlled  and  uncontrolled  ellipsoids.  Used  in  Subroutine 
COVER. 

CDET  =  Determinant  used  in  computing  volume  reduction  and  actual 
volume  of  ellipsoid  for  case  where  controls  are  used  and 
steady  state  is  assumed.  Used  in  Subroutine  COVER. 

CDEV1  =  Equivalent  to  (CBAR  -  VECMUC) ’  =  (C  -  uc)  '  . 

CDEV2  =  Equivalent  to  (CBAR  -  VECMUC)  =  (C  -  uc) . 

CEIGS  *  Complex  variable  counterpart  of  EIGS. 

CONST  =  Equivalent  to  (NUMREPS- 1) / (NUMREPS-MV- 1) ;  partial  value  used 
in  calculating  determinants  (DET)  in  m  ’best’  regressions 
mode.  Used  in  MAIN  program. 

CONST  1  =  Equivalent  to  (NUMREPS-2) / (NUMREPS (NUMREPS- 1 )) ;  partial  value 

used  in  computing  the  estimator  Sigma  Tilde  Hat.  Used  in 
Subroutine  COVER. 

C0NST2  =  Equivalent  to  ( I IN+ 1 )/ (NUMREPS (NUMREPS- 1 ))  ;  partial  value 
used  in  computing  the  estimator  Sigma  Tilde  Hat.  Used  in 
Subroutine  COVER. 

CONTROL  =  Character  vector  which  contains  names  of  controls  used  in  the 
evaluation . 

COVCV  *  Array  containing  covariance  matrix  of  controls. 

COVFILE  =  Name  of  datafile  containing  covariance  matrix,  if  it  is 
known . 

COVERAG  =  Array  containing  estimated  confidence  volume  coverage. 

COVERAG(l):  Controlled  coverage  on  steady  state  means, 
COVERAG ( 2) :  Uncontrolled  coverage  on  steady  state  means, 
COVERAG ( 3) :  Controlled  coverage  on  sample  mean,  and 
COVERAG (4) :  Uncontrolled  coverage  on  sample  mean. 

CTOP  =  Numerator  counterpart  of  CBOT. 

D1  =  Output  of  IMSL  Subroutine  LINV3F.  D1  is  one  of  two 

components  of  the  determininant  of  the  matrix  input  into 
LINV3F . 

D2  =  Output  of  IMSL  Subroutine  LINV3F.  D2  is  one  of  two 

components  of  the  determininant  of  the  matrix  input  into 
LINV3F . 

DET  =  Determinant  of  an  applicable  matrix.  Equivalent  to  D1»2*«D2. 

DEV  =  Equivalent  to  BETA*BUFF3  »  fl(C  -  ue>. 

DIFF  =  Equivalent  to  VC  -  VU. 

DUM  =  Dummy  array  for  use  in  IMSL  subroutine  calls,  when  output  of 

that  type  is  not  required.  Used  in  MAIN  program  and 
Subroutine  COVKNOW. 

DUMMY  =  Dummy  matrix  for  use  in  IMSL  subroutine  calls,  when  output  of 
that  type  is  not  required.  Used  in  Subroutine  COVER. 

EHAT  =  Matrix  containing  values  for  estimator  Sigma  Tilde  Hat. 

EIGS  =  Vector  variable  used  in  Subroutine  Cover  to  contain 
eigenvalues  derived  from  IMSL  subroutine. 

EPS  =  Covergence  criterion  used  as  input  to  IMSL  Subroutine  ZREAL2 

EXPL  =  Equivalent  to  BETA*BlTFF3  =  fi«Scy. 
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F 

F2 

FF 

FP 

FULL 

FULCOVC 

GAMM 

GAMMA 


HPH 

I 


11 

12 

I  AT 

IB 

IB2 

IBOT 
I  BUFF 


IC 

I COUNT 


I COVER 


ICTOT 


IER 

IH 


IH2 


II 
I  IN 


s  Used  in  computing  right  hand  side  of  confidence  region,  in 
Subroutine  COVER. 

*  Used  in  computing  right  hand  side  of  volume  reduction,  in 
Subroutine  COVER. 

a  Contains  F^  table. 

*  Equivalent  to  ROOT(l)««IP.  Used  in  Subroutine  FTABL  to  hold 
F-value  until  it  is  added  to  array  FF. 

=  Full  storage  mode  version  of  VCV. 

=  Full  covariance  matrix  of  all  controls  and  responses. 

=  Gamma  Hat  matrix,  used  in  Subroutine  COVER. 

=  Used  to  compute  actual  volume  of  controlled  ellipsoid,  in 
Subroutine  COVER. 

v  A/  <v 

=  Equivalent  to  T1»CDEV2  =  (C  -  uc)  -  uc)  . 

=  Counting  variable  used  in  DO  loops.  Also  used  to  count 
number  of  primary  variable  inputs  which  exceed  program 
parameters,  if  any. 

=  DO  loop  counting  variable,  used  in  Stepwise  section  of  MAIN 
program. 

=  DO  loop  counting  variable,  used  in  Stepwise  section  of  MAIN 
program. 

=  Tracks  which  regression  model  is  currently  being  evaluated. 

=  Index  of  source  block.  Used  in  Subroutine  GAUSS. 

=  Provides  same  function  as  IB.  Used  in  Stepwise  procedure  of 
MAIN  program. 

=  Denominator  for  values  computed  in  Functions  Cl  and  JCOMB. 

=  Array  containing  control  variables  in  ’best’  regression  model 
selected.  Identifies  variables  by  number,  not  model 
coef f icients . 

=  Counting  variable  used  in  DO  loops. 

=  Acts  as  reference  value,  in  Subroutine  COVER,  in  computing 
canonical  correlations  when  covariance  matrix  of  controls  is 
known. 

=  Indicator  array  of  coverage  for  a  particular  model;  0=No, 
l=Yes . 

ICOVER(l)  =  0,1  (Controls  present,  steady  state  assumed) 
IC0VER(2)  =  0,1  (No  controls,  steady  state  assumed) 

ICOVER (3)  =  0,1  (Controls  present,  Y(1000)) 

ICOVERU)  =0,1  (No  controls,  Y(  1000) ) 

=  Keeps  coverage  total  as  each  meta-experiment  is  performed. 
Used  in  MAIN  program  in  computing  average  coverage  over  all 
meta-experiments ;  COVERAG(I) = ICTOT (I) /META. 

=  Error  condition,  output  from  IMSL  subroutines  if  an  error 
condition  is  encountered. 

=  Inclusion  array,  of  both  controls  and  responses,  for  input  to 
IMSL  Subroutine  RLSUBM.  Used  in  MAIN  program  and  Subroutine 
COVER. 

=  Inclusion  array,  of  controls  only,  for  submatrix  of  selected 
model  used  as  input  to  IMSL  Subroutine  RLSUBM.  Used  in 
Subroutine  COVER. 

=  Counting  variable  used  in  DO  loops. 

=  Number  of  variables  in  current  model. 
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I  KNOW 


ILEFT 

I  LOOP 

IND1 

IND2 

INDEX 

INFILE 

I ORDER 

IP 

IQ 

IQQ 

IS 

IS2 

ITMAX 

I  TOP 
I  VAR 

IWRITE 

IX 

IZ 

J 

JJ 

JZ 

K 

KEEPERS 

KK 

KNX 

KP 

KZ 

L 

LAST 

LB 

M 

Ml 
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=  Flag  designating  whether  the  covariance  matrix  is  estimated 
(1KNOW=0) ,  or  known  (IKNOW=l). 

=  Number  of  pivots  left  to  process  in  computing  Function  C2; 
ILEFT  =  NFACT(IP)/(NFACT(IP-J)*NFACT(J)) . 

=  Equivalent  to  M-l  in  Function  NFACT.  Used  as  starting  point 
for  loop  which  computes  the  factorial  of  M. 

=  Array  for  keeping  ‘best’  evaluated  subset  of  size  J=1...NX. 
Used  in  Stepwise  procedure  of  MAIN  program. 

=  Array  for  maintaining  latest  subset  created  for  evaluation 
within  the  Stepwise  procedure  section  of  MAIN  orogram. 

=  Index  variable  for  SCBAR  and  SVECMU  when  computing  subvector 
of  the  control  means  in  Subroutine  COVER. 

=  Character  variable  used  in  initial  data  input  routine. 

INFILE  takes  name  of  input  file  if  datafile  option  is  chosen. 

=  Order  (number  of  rows)  of  input  or  output  matrix.  Used  in 
several  IMSL  Subroutines. 

=  Index  of  the  pivot  row  and  column,  used  in  Subroutine  GAUSS. 

=  Counting  variable  used  in  DO  loops. 

=  Counting  variable  used  in  DO  loops. 

=  Index  of  the  storage  block,  used  in  Subroutine  GAUSS. 

=  Provides  same  function  as  IS.  Used  in  Stepwise  procedure  of 
MAIN  program. 

=  Input  to  IMSL  Subroutine  ZREAL2,  defines  the  maximum  number 
of  iterations  to  use  in  finding  a  root. 

=  Numerator  for  values  computed  in  Functions  Cl  and  JCOMB. 

=  Integer  value  associated  with  the  individual  control 
variables  (i.e.  1=X1,  2=X2,  etc.). 

=  Flag  designating  whether  the  meta  experiment  mode  (IWRITE=0) 
or  best  ’m‘  regressions  mode  (IWRITE=1)  is  to  be  used. 

=  Equivalent  to  Z6 ,  used  as  input  to  IMSL  Subroutine  BECOVM. 

=  Counting  variable  used  in  DO  loops.  Used  in  initializing 
arrays  and  matrices  in  meta  loop  of  MAIN  program. 

=  Counting  variable  used  in  DO  loops. 

=  Counting  variable  used  in  DO  loops. 

=  Counting  variable  used  in  DO  loops.  Used  in  initializing 

arrays  and  matrices  in  meta  loop  of  MAIN  program. 

=  See  NX. 

=  Number  of  ’best’  regressions  to  keep  (See  M) . 

=  See  NUMREPS. 

=  Equivalent  to  2*»NX. 

=  Equivalent  to  k+1,  where  k  is  number  of  control  variates. 

Used  in  Subroutine  GAUSS. 

=  Counting  variable  used  in  DO  loops.  Used  in  initializing 
matrices  in  meta  loop  of  MAIN  program. 

=  Counting  variable  used  in  DO  loops. 

=  Variable  used  as  input  to  IMSL  Subroutine  ZREAL2,  contains 
initial  guess  of  root  for  defined  function. 

-  Equivalent  to  IP+1,  used  in  Subroutine  GAUSS. 

=  Counting  variable  used  in  DO  loops. 

=  Variable  used  as  input  to  IMSL  Subroutine  RLSUBM,  contains 
order  of  symmetric  matrix  stored  in  symmetric  mode. 


M2  =  Order  of  symmetric  matrix  SUBV.  Obtained  as  output  of  IMSL 

Subroutine  RLSUBM,  and  used  as  input  to  IMSL  Subroutine 
VCVTSF.  Used  in  Subroutine  COVER. 

MAX  =  Maximum  number  of  storage  matrices  in  Array  A,  for  storage  of 

matrices  created  by  calls  to  Subroutine  GAUSS.  Used  in 
Stepwise  section  of  MAIN  program. 

META  =  Number  of  meta  experiments.  This  is  required  input  data. 

METHOD  =  Flag  designating  whether  enumerated  subsets  (METH0D=0)  or 

stepwise  [forward  selection]  (METH0D=1)  method  is  to  be  used 
in  the  evaluation.  This  is  required  input  data. 

MM  =  Used  as  DO  loop  counting  variable  for  META  loop. 

MODELS  =  Saves  all  the  2N  enumerated  subset  models. 

MV  =  Counts  number  of  control  variables  contained  in  a  model  as 

defined  by  NK.  Used  in  MAIN  program  and  Subroutine  COVKNOW. 

N  =  Counting  variable  used  in  DO  loops. 

NBR  =  Vector  of  inputs  to  IMSL  Subroutine  BECOVM. 

NK  =  An  identification  array,  NK  is  a  binary  counter  with  a  list 

of  zeros  and  ones  which  indicate  the  presence  or  absence  of 
the  independent  variables  (i.e.  controls).  Also  note  that 
indexing  of  the  independent  variables  is  reversed. 

NROOT  =  Input  to  IMSL  Subroutine  ZREAL2,  defines  number  of  roots  to 
be  found. 

NSIG  *  Convergence  criterion  input  to  IMSL  Subroutine  ZREAL2.  A 

root  is  accepted  if  two  successive  approximations  to  a  given 
root  agree  in  the  first  NSIG  digits. 

NUMREG  =  Tracks  the  number  of  'best'  regressions,  when  that  mode  is 
used. 

NUMREPS  =  Number  of  replications  per  meta  experiment.  This  is  required 
input  data. 

NVAR  =  Total  number  of  variables  (i.e.  NVAR  =  NX  +  NY). 

NX  =  Number  of  candidate  control  variates.  This  is  required  data 

input . 

NY  =  Number  of  responses.  This  is  required  input  data. 

OBS  =  Equivalent  to  T2*YMD2. 

0BS2  =  Output  of  IMSL  Subroutine  VMULFF,  contains  product  of  the 

first  two  matrices  provided  as  input  to  the  subroutine.  Also 
used  in  determining  coverage  when  no  controls  are  used  and 
steady  state  is  assumed.  Used  in  Subroutine  COVER. 

OUTFILE  =  Character  variable  used  input  name  of  an  input  datafile 

created  during  manual  data  input,  if  desired.  Also  used  to 
input  name  for  file  to  contain  program  output. 

P  =  Output  probability,  of  IMSL  Subroutine  MDFD,  that  a  random 

variable  following  the  F-Distribution  with  degrees  of  freedom 
N1  and  N2  will  be  less  than  or  equal  to  input  Z.  Used  in 
Function  F.  This  value  supports  calculation  of  F-table  in 
Subroutine  FTABL. 

PI  =  Product  1,  numerator  used  in  computing  RNEXT,  in  Function  C2. 

PI  *  [PROD ( J=0 , IP)  (IQ+ (2* ( J+ 1 ) ) ) ]. 

P2  =  Product  2,  denominator  used  in  computing  RNEXT,  in  Function 

C2 .  P2  =  [PROD ( J=0 , IP)  (K-IQ- (2* J) ) ] . 

PI  =  Parameter  in  Subroutine  Cover  which  contains  value  of  pi. 
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P0VER2 

PROD 

RDET 

REGR 

REIGS 


RESPONS 

RHS 

RHS2 

RMIN 

RNEXT 

ROOT 


RSS 

SCBAR 


SIG 

SP 

SUBCOVC 

SUBV 

SUBVF 

SUM 

SUMDEV 


SUMVU 


SVECMU 

SYDOTC 

SYMCOVC 

T1 

T2 

TARGET 


Equivalent  to  NY/2.  Used  in  Subroutine  COVER  to  compute 
actual  volume  of  controlled  and  uncontrolled  ellipsoid. 

Keeps  cumulative  product  of  terms,  for  Functions  Cl,  C4,  and 
C5 . 

Array  of  determinants  of  matrices  associated  with  the 
specific  regression  models. 

Bookkeeping  array  for  best  M  regressions  to  keep.  For  array 
of  format  REGR(i,j,k);  A)  j  =  subset  size,  B)  k=l,  stores 
generalized  matrix;  and  k=2,  stores  pointer  to  model. 
Equivalent  to  Sqrt (EIGS (ICOUNT) ) ,  used  in  computing  the 
canonical  correlations  and  feasibility  bound  when  the 
covariance  matrix  of  controls  is  known.  Used  in  Subroutine 
COVER. 

Character  vector  used  to  contain  names  of  responses  used  in 
the  evaluation. 

Right  Hand  Side  of  the  confidence  region,  where  controls  are 
used,  per  Rao  (1967).  Used  in  Subroutine  COVER. 

Right  Hand  Side  of  the  confidence  region,  when  no  controls 
are  used  and  steady  state  assumed.  Used  in  Subroutine  COVER. 
Holds  minimum  Residual  Sums  of  Squares  (RSS)  values.  RSS 
values  used  in  determining  the  m  'best*  regressions,  when 
that  mode  is  selected. 

Equivalent  to  P1/P2  in  Function  C2.  Used  in  computing  TERM. 
Used  as  input/output  to  IMSL  Subroutine  ZREAL2,  in  Subroutine 
FTABL.  As  input,  contains  initial  guess  of  root.  As  output, 
contains  computed  root. 

Buffer  for  conditional  covariance  matrix. 

Subvector  of  CBAR,  used  to  find  vector  of  corrections  to 
control  Y  (variables  CDEV1  and  CDEV2) .  Used  in  Subroutine 
COVER. 

Level  of  Significance  associated  with  selection  criteria. 
Equivalent  to  RMIN,  used  in  write  statement  to  program  output 
f  ile . 

Submatrix  of  full  covariance  matrix  FULCOVC. 

Submodel  covariance  matrix  in. symmetric  storage. 

Full  storage  version  of  SUBV 
Keeps  sum  of  TERMs  in  Function  C2. 

Keeps  sum  of  DIFFs  for  each  meta-experiment. 

SUMDEV(l):  For  steady  state  means. 

SUMDEV(2):  For  sample  means. 

Keeps  sum  of  VUs  for  each  meta-experiment. 

SUMVU(l):  For  steady  state  means. 

SUMVU(2):  For  sample  means. 

Subvector  of  VECMUC,  used  to  find  vector  of  corrections  to 
control  Y.  Used  in  Subroutine  COVER. 

Equivalent  to  BUFF4-EXP  s  C2(Syy  -  fi«SCy) . 

Symmetric  storage  version  of  matrix  FULCOVC. 

W  mr 

Equivalent  to  CDEV1*SUBVF  =  (C  -  u^'tS,,"1. 

Equivalent  to  YMD1« SYDOTC . 

Matrix  used  in  computing  generalized  variance  of  Sigma  Tilde 
Hat  (i.e.  responses). 
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TEMP 


=  Input  vector  of  length  NBR(l),  used  in  IMSL  Subroutine 
BECOVM.  If  NBR(5)=0,  then  TEMP  must  contain  the  temporary 
means  when  NBR(4)=1.  Otherwise,  temp  is  work  storage. 

TERM  =  Generic  variable,  used  in  Functions  Cl  and  C2  to  save  results 
of  divisions  and  products,  respectively. 

TERM1  =  Partial  value  used  in  computing  volume  reduction  where 
controls  are  used  and  steady  state  assumed.  Used  in 
Subroutine  COVER. 

TERM2  =  Partial  value  used  in  computing  volume  reduction  where 
controls  are  used  and  steady  state  assumed.  Used  in 
Subroutine  COVER. 

TERM3  =  Partial  value  used  in  computing  volume  reduction  where 
controls  are  used  and  steady  state  assumed.  Used  in 
Subroutine  COVER. 

TIND  =  Array  for  keeping  current/temporary  'best'  variable  subset 
model  as  all  models  are  created  for  evaluation  in  Stepwise 
section  of  MAIN  program.  When  all  models  of  size  J=1,...,NX 
have  been  evaluated,  the  model  maintained  in  TIND  is  copied 
to  IND1,  as  the  'best*  subset  of  size  J. 

TITLE  =  Analysis  title,  input  during  initial  data  input  routine  and 
written  to  output  file. 

TMV  =  Tracks  number  of  variables  in  last  variable  subset  evaluated 

in  the  Stepwise  procedure  section  of  MAIN  program.  Acts  as  a 
flag  to  trigger  save  of  best  subset  of  a  certain  size. 

TOP  =  Numerator  portion  of  a  value  used  to  compute  C4,  in  Function 

C4 . 

TWO  «  Modified  version  of  matrix  determinant,  equivalent  to 

2«BIG«DET.  Used  as  determinant  bound  and  as  original  values 
for  elements  of  bookkeeping  array  REGR. 

UCDET  =  Determinant  used  in  computing  the  volume  reduction  where 

controls  are  used  and  steady  state  assumed,  and  the  actual 
volume  of  the  uncontrolled  ellipsoid.  Used  in  subroutine 
COVER. 

VC  *  Volume  of  controlled  ellipsoid,  used  in  Subroutine  COVER. 

VCV  =  Covariance  matrix  in  symmetric  storage  of  all  controls  and 

responses.  This  is  an  output  of  IMSL  Subroutine  BECOVM. 

VECCBAR  =  Vector  of  average  of  inputs  of  controls. 

VECMUC  =  Vector  of  theoretical  means  of  the  controls. 

VECMUY  =  Vector  of  steady  state  means  of  the  responses. 

VECYBAR  5  Vector  of  sample  means  of  the  responses. 

VOLRED  ■  VOLRED(l)  is  volume  reduction  due  to  controls;  V0LRED(2)  is 
not  used. 

VR  =  Array  of  Variance  Reduction  values,  derived  from  calulating 

selection  criterion  for  each  regression  model. 

VU  =  Volume  of  uncontrolled  ellipsoid,  used  in  Subroutine  COVER. 

WK  =  Array  used  by  IMSL  Subroutines  EIGRF.  Provides  work  area  for 

subroutine  to  use  in  performing  its  function. 

WXAREA  =  Array  used  by  IMSL  Subroutine  LINV3F.  Provides  same  function 
as  WX. 

X  =  Data  matrix  for  a  single  meta-experiment. 

XFILE  *  Datafile  containing  [controls ! response ]  data. 
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XM 

YBAR 

YBHAT 

YMD1 

YMD2 

YMD3 

YMD4 

Z 

Z1 

Z2 

Z3 

Z4 

Z5 

Z6 

Z7 

Z8 

NOTES:  ( 


=  Output  vector,  of  length  NBR(l),  of  IMSL  Subroutine  BECOVM. 
This  vector  contains  the  variable  means. 

=  Sample  mean  vector. for  controls.. 

=  Equivalent  to  Y  -  fi(C  -  uc)  =  Y(£). 

=  Equivalent  to  (Y(fi)  -  Uy) ’ . 

=  Equivalent  to  (Y(fi)  -  Uy) ' . 

=  Equivalent  to  (YBAR  -  VECMUY) 1  =  (Y  -  uy)  ’ . 

=  Equivalent  to  (YBAR  -  VECMUY)  =  (Y  -  uJ . 

=  Input  constant,  to  IMSL  Subroutine  MDFD,  to  which  integration 
is  performed.  Z  must  be  greater  than  or  equal  to  zero.  Used 
in  Function  F. 

*  Program  parameter  defining  maximum  number  of  Control 
Variables  (See  NX)  the  program  is  set  to  handle.  This 
parameter  is  also  used  in  Subroutine  Cover. 

=  Program  parameter  defining  maximum  number  of  Control 
Responses  (See  NY)  the  program  is  set  to  handle.  This 
parameter  is  also  used  in  Subroutine  Cover. 

=  Program  parameter  equivalent  to  Z1+Z2.  This  parameter  is 
also  used  in  Subroutine  Cover. 

=  Program  parameter  defining  maximum  number  of  'best' 
regressions  which  may  be  kept. 

=  Program  parameter  equivalent  to  2»*Z1.  This  parameter  also 
used  in  Subroutine  Cover. 

=  Program  parameter  defining  maximum  number  of  replications  per 
meta-experiment  allowed. 

=  Program  parameter  defining  maximum  number  of  meta  experiments 
allowed. 

=  Program  parameter  equivalent  to  (Z3* (Z3+1) ) /2.  This 
parameter  also  used  in  Subroutine  Cover. 


1)  S 


/  \ 

S'  Q  ' 

cc  1  °cy  1 


Syc  •  Syy 
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VARIABLE  PARAMETERS  MATRIX 
(Continued) 


VARIABLE 

I  NT 

REAL 

fariab] 

CHAR 

Le  Typt 
COM 

LOGIC 

ARRAY 

Precj 

S 

sion 

D 

Size/Comments 

CC1 

X 

X 

CC2 

X 

X 

CC3 

X 

X 

CC4 

X 

X 

CDET 

X 

X 

CDEV1 

X 

X 

X 

lxZl 

CDEV2 

X 

X 

X 

Zlxl 

CEIGS 

X 

X 

Z2 

CONST 

X 

■■ 

X 

CONST 1 

im 

X 

B 

X 

C0NST2 

X 

BB 

X 

CONTROL 

X 

X 

C : 25  A:Z1 

COVC  V 

BB 

X 

X 

X 

ZlxZl 

COVFILE 

II 

X 

25 

COVERAG 

X 

X 

X 

4 

CTOP 

X 

X 

D1 

X 

X 

D2 

X 

X 

DET 

X 

X 

DEV 

X 

X 

X 

Z2xl 

DIFF 

X 

X 

DUM 

X 

X 

X 

Z2 

DUMMY 

X 

X 

X 

Z2xZ2 
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VARIABLE  PARAMETERS  MATRIX 
(Continued) 


VARIABLE  !  I NT  1 REAL 


Variable  Type  ! 

CHAR  !  COM  ! LOGIC! ARRAY! 


Precision 
S  !  D 


! Size/Comments ! 


EHAT 

X 

X 

EIGS 

X 

X 

EPS 

X 

EXPL 

X 

X 

F 

X 

F2 

X 

FF 

X 

X 

FP 

X 

FULL 

X 

X 

FULCOVC 

X 

X 

GAMM 

X 

X 

HPH 

X 

X 

I 

X 

11 

X 

12 

X 

I  AT 

X 

IB 

X 

IB2 

X 

IBOT 

X 

I  BUFF 

X 

X 

IC 

X 

I COUNT 

X 

I COVER 

;  X 
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VARIABLE  PARAMETERS  MATRIX 
(Continued) 


VARIABLE 

I  NT 

REAL 

/ariab 

CHAR 

ie  Typ< 
COM 

a 

LOGIC 

ARRAY 

Preci 

S 

ision 

D 

Size/Comments 

ICTOT 

X 

X 

4 

IER 

X 

IH 

X 

X 

Z3 

IH2 

X 

X 

Z1 

II 

X 

I  IN 

X 

I  KNOW 

X 

ILEFT 

X 

I  LOOP 

X 

IND1 

X 

X 

Zl  +  1 

IND2 

X 

X 

Zl  +  1 

INDEX 

X 

INFILE 

X 

25 

I ORDER 

X 

IP 

X 

IQ 

X 

IQQ 

X 

IS 

X 

IS2 

X 

ITMAX 

X 

I  TOP 

X 

I  VAR 

X 

IWRITE 

X 
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VARIABLE  PARAMETERS  MATRIX 
(Continued) 


VARIABLE 

I  NT 

1 

REAL 

fariabl 

CHAR 

e  Typ« 
COM 

LOGIC 

ARRAY 

Preci 

S 

M  —  . 

sion 

D 

Size/Comments 

IX 

X 

IZ 

X 

J 

X 

_ _ 

JJ 

X 

JZ 

X 

K 

X 

KEEPERS 

X 

KK 

X 

_  _ 

KNX 

X 

KP 

X 

KZ 

X 

L 

X 

LAST 

m 

X 

X 

LB 

■■ 

X 

M 

X 

Ml 

X 

M2 

X 

MAX 

X 

Program 

Parameter 

META 

X 

METHOD 

X 

MM 

X 

MODELS 

X 

X 

Z5xZl 
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VARIABLE  PARAMETERS  MATRIX 
(Continued) 


VARIABLE 

I  NT 

REAL 

Jaridbi 

CHAP 

le  Typ« 
COM 

_ 

LOGIC 

ARRAY 

Precj 

S 

LSion 

D 

Size/Comments 

MV 

X 

1  | 

N 

X 

NBR 

X 

X 

6 

NK 

X 

X 

Z3 

NROOT 

X 

NSIG 

X 

NUMREG 

X 

NUMREPS 

X 

NVAR 

X 

NX 

X 

NY 

X 

OBS 

X 

X 

X 

lxl 

0BS2 

X 

X 

X 

lxl 

OUTFILE 

X 

25 

P 

X 

X 

PI 

X 

X 

P2 

X 

X 

PI 

X 

X 

PROD 

X 

X 

RDET 

X 

X 

REGR 

X 

X 

X 

Z4xZlx2 

REIGS 

X 

X 

X 

Z2 

RESPONS 

X 

X 

C: 25  A:Z2 
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VARIABLE  PARAMETERS  MATRIX 
(Continued) 


VARIABLE 

I  NT 

\ 

REAL 

!  ariabi 
CHAR 

e  Typt 
COM 

LOGIC 

ARRAY 

Precj 

S 

Sion 

D 

Size/Comments 

RHS 

X 

X 

RHS2 

X 

X 

RMIN 

X 

X 

RNEXT 

X 

X 

ROOT 

X 

X 

X 

1 

RSS 

X 

X 

X 

Z2xZ2 

SCBAR 

X 

X 

X 

Z1 

SIG 

X 

X 

SP 

X 

X 

SUBCOVC 

X 

X 

X 

( Z 1 (Zl+1) ) /2 

SUBV 

X 

X 

X 

Z8 

SUBVF 

X 

X 

X 

Z3xZ3 

SUM 

V 

X 

SUMDEV 

X 

X 

X 

2 

SUMVU 

X 

X 

X 

2 

SVECMU 

X 

X 

X 

Z1 

SYDOTC 

X 

X 

X 

Z2xZ2 

SYMCOVC 

X 

X 

X 

( Z 1 (Z 1 + 1 ) ) /2 

T1 

X 

X 

X 

lxZl 

T2 

X 

X 

X 

lxZ2 

TARGET 

X 

X 

X 

Z2xZ2 

TEMP 

X 

X 

X 

Z3 

TERM 

X 

X 

1 1 1 


VARIABLE  PARAMETERS  MATRIX 
(Continued) 


!  Variable  Type  !  Precision  !  ! 

VAP.  .ib  >E  !  INT  IREAL  ICHAR  !  COM  ! LCGIC !  ARRAY S  !  D  !  Size/Comments ! 


VARIABLE  PARAMETERS  MATRIX 
(Continued) 


Variable  Type 


Precision 


VARIABLE 

I  NT 
_ _ _ 

REAL 

CHAR 

COM 

LOGIC 

ARRAY 

s 

D 

Size/Comments 

YBHAT 

■■ 

X 

X 

X 

Z2 

fMDl 

H 

X 

X 

X 

lxZ2 

YMD2 

BB 

X 

X 

X 

Z2xl 

YMD3 

X 

X 

X 

lxZ2 

YMD4 

X 

X 

X 

Z2xl 

Z 

X 

X 

Z1 

X 

Program 

Parameter 

Z2 

X 

Program 

Parameter 

Z3 

X 

Program 

Parameter 

Z4 

X 

Program 

Parameter 

Z5 

X 

Program 

Parameter 

Z6 

X 

Program 

Parameter 

Z7 

X 

Program 

Parameter 

Z8 

X 

Program 

Parameter 

_ _ _ 

_ _ 

_ _ 
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V.  SUBROUTINE  LISTING  AND  DESCRIPTION 


BECOVM  =  An  IMSL  subroutine  which  computes  means  and  variance- 
covariance  matrix. 

COVER  =  This  subroutine  does  the  coverage  and  volume  reduction 
calculations  for  the  optimal  control  subset. 

COVKNOW  =  This  subroutine  returns  the  generalized  variance  of  sigma 
tilde  hat. 

EIGRF  =  An  IMSL  subroutine  which  computes  eigenvalues  and 

(optionally)  eigenvectors  of  a  real  general  matrix  in  full 
storage  mode. 

FTABL  =  This  subroutine  computes  an  F  table,  to  the  power  P. 

GAMMA  =  An  implicit  FORTRAN  function  which  provided  Gamma  Function 

value . 

GAUSS  =  This  subroutine  performs  the  pivots  for  variable 

introduction  into  regression  models.  (Furnival  and  Wilson 
1974) 

KEEPIT  =  This  subroutine  finds  the  model  of  a  candidate  regression. 

LINV3F  =  An  IMSL  subroutine  which  computes  in-place  inverse, 

equation  solution,  and/or  determinant  evaluation,  uses 
full  storage  mode. 

MDFD  =  An  IMSL  subroutine  which  computes  the  F  probability 

distribution  function. 

RLSUBM  =  An  IMSL  subroutine  which  performs  retrieval  of  a  symmetric 
submatrix  from  a  matrix  stored  in  symmetric  mode  by  RLSTP. 

VCVTFS  =  An  IMSL  subroutine  which  performs  storage  mode  conversion 
of  matrixes  (full  to  symmetric) . 

VCVTSF  -  An  IMSL  subroutine  which  performs  storage  mode  conversion 
of  matrixes  (symmetric  to  full). 

VMULFF  =  An  IMSL  subroutine  which  performs  matrix  multiplication 
(full  storage  mode) . 

ZREAL2  =  An  IMSL  subroutine  which  computes  the  real  zeros  of  a  real 
function  -  to  be  used  when  initial  guesses  are  good. 


Note:  The  following  subroutine  is  not  used  in  program,  but  is 
referenced  in  description  of  subroutine  RLSUBM.  In  addition,  complete 
IMSL  subroutine  usage  descriptions  are  provided  in  Attachment  D. 


[RLSTP ] 


=  An  IMSL  subroutine  which  performs  regression  model 

selection  using  a  forward  stepwise  algorithim,  with  results 
available  after  each  step. 
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Cl 


C2 


C3 


C4 


C5 


CFRONT 


F 


JCOMB 


NFACT 


VI.  FUNCTION  LISTING  AND  DESCRIPTION 


=  Equivalent  to  [PROD (1=1 , IP)  ( (K-IQ-I) / (K (K-IQ-1) ) ]. 

Partial  value  used  in  computing  selection  criterion  for  a 
model,  when  the  covariance  matrix  of  controls  is  estimated. 
Used  in  Function  C3. 

=  Equivalent  to  (1  ♦  [SUM(J=1,P)  JCOMB (P . J) • (Q(Q+2) ... (Q+ 

2 (J- 1) ) )  /  ( (K-Q-2) . . . (K-Q-2J) ) ]} .  Partial  value  used  in 
computing  the  selection  criterion  when  the  covariance 
matrix  is  estimated.  Used  in  Function  C3. 

=  Equivalent  to  Cl  •  C2.  Partial  value  used  in  computing 
selection  criterion  for  a  model,  when  covariance  matrix  of 
controls  is  estimated.  This  provides  a  loss  factor  for  the 
measure  of  efficiency  of  control  variables  as  defined  by 
Rubenstein  and  Marcus  (1985).  Used  in  MAIN  program. 

=  Equivalent  to  [PR0D(I=1 , IP)  ( (K-IQ-I) / (K-IQ-1) )] .  Partial 
value  used  in  computing  selection  criterion  for  a  model, 
when  covariance  matrix  of  controls  is  known.  Used  in  MAIN 
program. 

=  Equivalent  to  [PR0D(I=1 ,IP)  ( (KK-IQ-I) / (KK-IQ-1) ) ] . 

Partial  value  used  in  computing  selection  criterion  for  a 
model,  regardless  of  whether  covariance  matrix  is  known  or 
estimated.  Used  in  MAIN  program. 

=  Equivalent  to  I«[ (K-IQ-1) / (K-IQ-IP) ] .  This  value  is  used 
in  calculating  the  100(l-a)X  confidence  ellipsoid  about  the 
responses  as  defined  by  Rao  (1967).  Used  in  MAIN  program. 

=  A  single-arguement  real  function  subprogram  used  by  IMSL 
Subroutine  ZREAL2.  F  defines  the  function  for  which  the 
roots  are  to  be  found. 

=  Equivalent  to  [NFACT(N) / (NFACT(N-M) *NFACT(M) ) ] ,  which  is 
the  number  of  possible  combinations  of  N  items  taken  M  at  a 
time . 

=  Computes  factorial  of  X  (i.e.  X!  *  X«(X-l)*(X-2)«. . .«1) . 
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ATTACHMENT  A:  VARIABLE  SUBSET  SELECTION  PROGRAM  CODE  (FORTRAN) 

This  program  listing  is  already  provided  as  Appendix  A  to  the 
thesis . 
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ATTACHMENT  B:  FORMAT  FOR  DATA  INPUT  FILES 


In  the  event  that  program  data  will  be  input  through  the  use  of  a 
datafile,  the  following  format  is  required.  It  should  be  noted  that  the 
data  can  be  placed  as  shown  or  on  seperate  lines  for  each  value  see 
examples  1  and  2,  respectively),  the  order  remains  as  shown.  In 
addition,  a  combination  of  these  two  formats  may  be  used,  if  desired. 


NX  NY  KEEPERS  NUMREPS  META 

I KNOW  I WRITE  METHOD 

SIG 

VECMUC  (NX) 

VECCBAR  (NX) 

VECMUY  (NY) 

VECYBAR  (NY) 

TITLE 

CONTROL  (NX) 

RESPONS  (NY) 

CCOVFILE;  IF  IKN0W=1> 

<COVCV  (NXxNX) ;  IF  COVFILE=INFILE  AND  IKN0W=1> 
XFILE 


NOTES : 

U]  VARIABLE  NAME  *  SPECIFIC  DATA  TO  PLACE  AT  THAT  POINT. 

[2]  (XI)  =  XI  ROWS,  OR  TOTAL  OF  XI  ELEMENTS,  EXPECTED. 

[3]  (XlxX2)  =  XI  ROWS  BY  X2  COLUMNS,  OR  TOTAL  OF  XlxX2  ELEMENTS, 
EXPECTED. 

[4]  <X;Y>  =  COMMENTS  ON  OPTIONAL  FILE  DATA.  X  IS  REQUIRED  DATA  IF 
CONDITION  Y  IS  MET,  OTHERWISE  DO  NOT  INSERT  THIS  DATA  INTO  THE 
DATAFILE. 

[5]  DO  NOT  PLACE  COMMAS  BETWEEN  DATA  VALUES. 


DATAFILE  FORMAT  EXAMPLES: 

Example  1_: 

5  2  1  10  100 

0  0  0 

0.90 

0  0  0  0  0 

0  0  0  0  0 

0  0 

0  0 

EXAMPLE. 1 

Cl  C2  C3  C4  C5 
R1  B2 

SUMMARY  EX*1 


lie 


5 

2 

1 

10 

100 

0 

0 

0 

0.90 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

EXAMPLE  2 

Cl 

C2 

C3 

C4 

C5 

HI 

R2 

SUMMARY  EX* 2 


ATTACHMENT  C:  PROGRAM  FLOW  DIAGRAMS 
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Program  Elements  Called  By  Subroutine  COUER 


FTflBL 
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Program  Elements  Called  By  COUKNOM,  FTflBL,  and  C3 


ATTACHMENT  D 


IMSL  SUBROUTINE  USAGE  DESCRIPTIONS 
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IMSL  ROUTINE  NAME  -  BECOVM 

PURPOSE  -  Merns  and  Variance-Covariance  matrix 

USAGE  -  CALL  BECOVM (X , IX , NBR, TEMP , XM, VCV , IER) 

ARGUEMENTS  X  -  ON  INPUT:  X  is  an  NBR(3)  by  NBR(l)  suboatrix  of  the 

matrix  (call  it  XX)  of  data  for  which  means, 
variances  and  covariances,  or  corrected  sums  of 
squares  and  cross  products  are  desired.  The 
last  submatrix  in  XX  may  have  fewer  than  NBR(3) 
rows . 

ON  OUTPUT:  The  rows  of  X  have  been  adjusted  by  the 
temporary  means. 

IX  -  Input,  row  dimension  of  X  exactly  as  specified  in 
the  dimension  statement  in  the  calling  program. 

NBR  -  Input  vector  of  length  6.  NBR(I)  contains,  when 
1=1,  Number  of  variables. 

1=2,  Number  of  observations  per  variable  in  XX. 
1=3,  Number  of  observations  per  variable  in  each 
submatrix  X,  not  including  the  last  submatrix 
where  the  number  may  be  less  than  or  equal  to 
NBR(3).  However,  NBR(3)  should  be  the  same 
for  all  calls. 

1=4,  The  number  of  the  submatrix  stored  in  X. 
1=5,  The  temporary  mean  indicator.  If  NBR(5)=0, 
the  user  supplies  temporary  means  in  TEMP. 
Otherwise,  the  first  row  of  XX  (or  first  row 
of  X  when  NBR(4)=1)  is  utilized. 

1=6,  The  VCV  option.  If  NBR(6)=0,  VCV  contains 
the  Variance-Covariance  matrix.  Otherwise, 
VCV  contains  the  corrected  cross  sums  of 
squares  and  cross-products  matrix. 

TEMP  -  Input  vector  of  length  NBR ( 1 ) .  If  NBR(5)=0,  TEMP 
must  contain  the  temporary  means  when  NBR(4)=1. 
Otherwise,  TEMP  is  work  storage. 

XM  -  Output  vector  of  length  NBR(l)  containing  the 
variable  means. 

VCV  -  Output  NBR(l)  by  NBR(l)  matrix  stored  in  symmetric 
storage  mode  requiring  (NBR( 1) » (NBR( 1) +1) ) /2 
storage  locations.  VCV  contains  the  Variance- 
Covariance  matrix  or  the  corrected  sums  of 
squares  and  cross-products  matrix,  as  controlled 
by  the  VCV  option,  NBR (6) . 
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IER  -  Error  parameter.  (Output) 

Terminal  Error: 

IER=129,  Indicates  that  N3R(4)  is  less  than  1  or 
NBR(3) « (NBR(4) -1)  exceeds  NBR(2) . 

I ER= 130,  Indicates  that  NBR(l)  is  less  than  1  or 
NBR(2)  is  less  than  2  or  NBR(3)  exceeds 
NBR(2) . 

PRECISION/HARDWARE  -  Single  and  Double/H32 

-  Single/H36 , H48 , H60 

REQD  IMSL  ROUTINES  -  UERTST ,  UGETIO 

NOTATION  -  Information  on  special  notation  and  conventions  is 

available  in  the  manual  or  through  IMSL  routine 
UHELP . 
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IMSL  ROUTINE  NAME 


EIGRF 


PURPOSE 

USAGE 

ARGUEMENTS 


-  Eigenvalues  and  (optionally)  eigenvectors  of  a  real 
general  matrix  in  full  storage  mode. 

-  CALL  EIGRF(A.N,IA,IJOB,W,Z,IZ,WK,IER) 


A  -  The  input  real  general  matrix  of  order  N  whose 

eigenvalues  and  eigenvectors  are  to  b  ecomputed. 
Input  A  is  destroyed  if  IJOB  is  equal  to  0  or  1. 

N  -  The  input  order  of  the  matrix  A. 

IA  -  The  input  row  dimension  of  matrix  A  exactly  as 

specified  in  the  dimension  statement  in  the 
calling  program. 

IJOB  -  The  input  option  parameter.  When 

IJOB=0,  Compute  eigenvalues  only. 

IJ0B=1,  Compute  eigenvalues  and  eigenvectors. 
IJ0B=2,  Compute  eigenvalues,  eigenvectors,  and 
performance  index. 

IJ0B=3,  Compute  performance  index  only.  If  the 
performance  index  is  computed,  it  is  returned 
in  WK(1).  The  routines  have  performed  (Well, 
Satisfactorily,  Poorly)  if  WK(1)  is  (Less 
than  1,  Between  1  and  100,  greater  than  100). 

W  -  The  output  complex  vector  of  length  N,  containing 

the  eigenvalues  of  A. 

Z  -  The  output  N  by  N  complex  matrix  containing  the 

eigenvectors  of  A.  The  eigenvector  in  column  J 
of  Z  corresponds  to  the  eigenvalue  W(J) .  If 
IJ0B=0,  Z  is  not  used. 


IZ 


WK 


-  The  input  row  dimension  of  matrix  Z  exactly  as 
specified  in  the  dimension  statement  in  the 
calling  program.  IZ  must  be  greater  than  or 
equal  to  N  if  IJOB  is  not  equal  to  zero. 


Work  area,  the  length  of  WK  depends  on  the  value  of 
IJOB,  when 


IJ0B=0 , 

The 

length 

of 

WK 

is 

at 

least 

N. 

I  J0B=  1 , 

The 

length 

of 

WK 

is 

at 

least 

2N. 

I J0B=2 , 

The 

length 

of 

WK 

is 

at 

least 

(2+N) N 

IJ0B=3 , 

The 

length 

of 

Wk 

is 

at 

least 

1. 
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IER  -  Error  parameter.  (Output) 

Terminal  Error: 

IER=128+J,  Indicates  that  EQRH3F  failed  to 
converge  on  eigenvalue  J.  Eigenvalues  J+l, 

J  +  2,  ....  N  have  been  computed  correctly. 
Eigenvalues  1,...,J  are  set  to  zero. 

IER=1  or  2,  Eigenvectors  are  set  to  zero.  The 
performance  index  is  set  to  1000. 

Warning  Error  (with  fix) 

IER=66,  Indicates  IJOB  is  less  than  0  or  IJOB  is 
greater  than  3.  IJOB  set  to  1. 

IER=67,  Indicates  IJOB  is  not  equal  to  zero,  and 
IZ  is  less  than  the  order  of  matrix  A.  IJOB 
is  set  to  zero. 

PRECISION/HARDWARE  -  Single  and  Double/H32 

-  Single/H36 ,H48 , H60 

REQD  IMSL  ROUTINES  -  EBALAF ,  EBBCKF ,  EHBCKF ,  EHESSF ,  EQRH3F ,  UERTST , 

UGETIO 

NOTATION  -  Information  on  special  notation  and  conventions  is 

available  in  the  manual  introduction  or  through 
IMSL  routine  UHELP. 
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IMSL  ROUTINE  NAME  -  LINV3F 

PURPOSE  -  In  place  inverse,  equation  solution,  and/or 

determinant  evaluation  -  full  storage  mode. 

USAGE  -  CALL  LINV3F (A, B , I JOB ,N, IA.D1 ,D2 , WXAREA , IER) 

ARGUEMENTS  A  -  Input/output  matrix  of  dimension  N  by  N.  See 

parameter  1J0B. 

B  -  Input/output  vector  of  length  N  when  IJ0B=2  or  3. 

Otherwise,  B  is  not  used. 

On  input:  B  contains  the  right  hand  side  of  the 
equation  AX  =  B. 

On  output:  The  solution  X  replaces  B. 

IJOB  -  Input  option  parameter.  IJ0B=I  implies: 

1=1,  Invert  matrix  A.  A  is  replaced  by  its 
inverse . 

1=2,  Solve  the  equation  AX=B.  A  is  replaced  by 
the  LU  decomposition  of  a  rowwise  permutation 
of  A,  where  U  is  upper  triangular  and  L  is 
lower  triangular  with  unit  diagonal.  The 
unit  diagonal  of  L  is  not  stored. 

1=3,  Solve  AX=B  and  invert  matrix  A. 

1=4,  Compute  the  determinant  of  A.  A  is 

replaced  by  the  LU  decomposition  of  a  rowwise 
permutation  of  A. 

N  -  Order  of  A.  (Input) 

IA  -  Row  dimension  of  matrix  A  exactly  as  specified  in 
the  dimension  statement  in  the  calling  program. 
(Input) 

D1  -  Input/Output.  If  the  D1  and  D2  components  of 
D2  determinant(A)  =  D1*2((D2  are  desired,  input 

D1.GE.0.  Otherwise,  input  D1.LT.0.  D2  is  never 
input. 

WKAREA  -  Work  area  of  length  at  least  2«N  for  IJOB=l  or  3. 

Work  area  of  length  at  least  N  for  IJ0B=2  or  4. 

IER  -  Error  parameter.  (Output) 

Terminal  Error: 

IER=130,  Indicates  that  matrix  A  is 
algorithmically  singular. 

Warning  With  Fix: 

IER=65,  Indicates  that  IJOB  was  less  than  1  or 
greater  than  4.  IJOB  is  assumed  to  be  4. 
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PRECISION/HARDWARE  -  Single  and  Double/H32 

-  Single/H36 ,H48 ,H60 

REQD  IMSL  ROUTINES  -  LUDATN ,  LUELMN,  UERTST ,  UGETIO 

NOTATION  -  Information  on  special  notation  and  conventions  is 

available  in  the  manual  introduction  or  through 
IMSL  routine  UHELP . 
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IMSL  SOUTINE  NAME 

PUBPOSE 

USAGE 

ASGUEMENTS  F 

N1 

N2 

P 

IER 


PRECISION /HARDWARE 
REQD  IMSL  ROUTINES 
NOTATION 


MDFD 

F  probability  distribution  function 
CALL  MDFD(F,N1 ,N2,P,IER) 

Input  constant  to  which  integration  is  performed. 

F  must  be  greater  than  or  equal  to  zero. 

Input  first  degree  of  freedom.  A  positive  integer. 

Input  second  degree  of  freedom.  A  positive 
integer . 

Output  probability  that  a  random  variable  following 
the  F  distribution  with  degrees  of  freedom  N1  and 
N2  will  be  less  than  or  equal  to  input  F. 

Error  parameter.  (Output) 

Terminal  Error: 

IER=129,  Indicates  either  N1  or  N2  is  less  than 
one  or  N2+N2  is  greater  than  20,000.  P  is 
set  to  positive  machine  infinity. 

IER=130,  Indicates  F  is  less  than  zero.  P  is 
set  to  positive  machine  infinity. 

Single/All 

MERRC=ERFC ,  UERTST ,  UGETIO 

Information  on  special  notation  and  conventions  is 
available  in  the  manual  introduction  or  through 
IMSL  routine  UHELP. 
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IMSL  ROUTINE  NAME  -  RLSUBM 


PURPOSE 

USAGE 

ARGUEMENTS  A 

M 

IH 


S 

N 

PRECISION/ HARD WARE 

REQD  IMSL  ROUTINES 
NOTATION 


Retrieval  of  a  symmetric  submatrix  from  a  matrix 
stored  in  symmetric  storage  mode  by  RLSTP. 

CALL  RLSUBM(A.M,IH,S,N) 

M  by  M  symmetric  matrix  stored  in  symmetric  mode. 
(Input)  A  is  a  vector  of  length  M«(M+l)/2. 

Order  of  the  matrix  A.  (Input) 

Vector  of  length  M.  (Input) 

If  IH(I) =IH(J)  =  1  where  J  and  1  =  1,2 . M,  the 

(I,J)-th  element  of  A  will  be  included  in  the 
submatrix  S. 

Otherwise,  the  (I,J)-th  element  of  A  will  not  be  in 
the  submatrix  S  on  output. 

Symmetric  submatrix  of  matrix  A.  (Output) 

S  is  a  vector  of  length  N»(N+l)/2.  A  and  S  may 
share  the  same  storage  if  it  is  not  necessary  to 
retain  the  original  matrix.  See  remarks  for 
RLSTP. 

Order  of  the  submatrix  S.  (Output) 

Single  and  Double/H32 
Sing 1 e/H36 , H48 , H60 

None  required. 

Information  on  special  notation  and  conventions  is 
available  in  the  manual  introduction  or  through 
IMSL  routine  UHELP. 
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IMSL  ROUTINE  NAME  -  VCVTFS 

PURPOSE  -  Storage  mode  conversion  of  matrices  (Full  to 

Symmetric) 

USAGE  -  CALL  VCVTFS (A ,N , IA,B) 

ARGUEMENTS  A  -  Input  matrix  of  dimension  N  by  N.  A  contains  a 

symmetric  matrix  stored  in  full  mode. 

N  -  Order  of  matrix  A.  (Input) 

IA  -  Row  dimension  of  matrix  A  exactly  as  specified  in 
the  dimension  statement  in  the  calling  program. 
(Input) 

B  -  Output  vector  of  dimension  NMN+D/2  containing 

matrix  A  in  symmetric  storage  mode. 

PRECISION /HARDWARE  -  Single  and  Double/H32 

-  Single/H36,H48,H60 

REQD  IMSL  ROUTINES  -  None  required. 

NOTATION  -  Information  on  special  notation  and  conventions  is 

available  in  the  manual  introduction  or  through 
IMSL  routine  UHELP . 
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IMSL  ROUTINE  NAME 


VCVTSF 


PURPOSE 

USAGE 

ARGUEMENTS  A 

N 

B 

IB 

PRECISION/HARDWARE 

REQD  IMSL  ROUTINES 
NOTATION 


Storage  mode  conversion  of  matrices  (Symmetric  to 

Full) 

CALL  VCVTSF (A,N,B , IB) 

Input  vector  of  length  N*(N+l)/2  containing  an  N  by 
N  symmetric  matric  stored  in  symmetric  storage 
mode . 

Order  of  matrix  A.  (Input) 

Output  matrix  of  dimension  N  by  N  containing  matrix 
A  in  full  storage  mode. 

Row  dimension  of  matrix  B  exactly  as  specified  in 
the  dimension  statement  in  the  calling  program. 
(Input) 

Single  and  Double/H32 

Single/H36 , H48 , H60 

None  required. 

Information  on  special  notation  and  conventions  is 

available  in  the  manual  Introduction  or  through 

IMSL  routine  UHELP. 
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IMSL  SOUTINE  NAME 

PURPOSE 

USAGE 

ARGUEMENTS  A 
B 
L 
M 

N 

IA 

IB 

C 

IC 

IER 

PRECIS ION/ HARDWARE 

REQD  IMSL  ROUTINES 
NOTATION 


VMULFF 

Matrix  multiplication  (Full  storage  mode) 

CALL  VMULFFU.B.L.M.N.IA.IB.C.IC.IEB) 

L  by  M  matrix  stored  in  full  storage  mode.  (Input) 
M  by  N  matrix  stored  in  full  storage  mode.  (Input) 
Number  of  rows  in  A.  (Input) 

Number  of  columns  in  A  (same  as  number  of  rows  in 
B) .  (Input) 

Number  of  column  in  B.  (Input) 

Row  dimension  of  matrix  A  exactly  as  specified  in 
the  dimension  statement  in  the  calling  program. 
(Input) 

Bow  dimension  of  matrix  B  exactly  as  specified  in 
the  dimension  statement  in  the  calling  program. 
(Input) 

L  by  N  matrix  containing  the  product  C  =  A*B. 
(Output) 

Row  dimension  of  matrix  C  exactly  as  specified  in 
the  dimension  statement  in  the  calling  program. 
(Input) 

Error  parameter.  (Output) 

Terminal  Error: 

IER=129,  Indicates  A,  B,  or  C  was  dimensioned 
incorrectly. 

Single  and  Double/H32 
Single/H36 ,H48 ,H60 

UERTST ,  UGETIO 

Information  on  special  notation  and  conventions  is 
available  in  the  manual  introduction  or  through 
IMSL  routine  UHELP. 
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IMSL  ROUTINE  NAME  -  ZREAL2 

PURPOSE  -  Determine  the  real  zeros  of  a  real  function  -  to  be 

used  when  initial  guesses  are  good. 

USAGE  -  CALL  ZREAL2 (F ,EPS , EPS2 , ETA , NSIG ,N , X , ITMAX , IER) 

ARGUEMENTS  F  -  A  single-arguement  real  function  subprogram 

supplied  by  the  user.  (Input) 

EPS  -  Convergence  criterion.  (Input) 

A  root,  X(I),  is  accepted  if  ABS(X(I) ) .LE.EPS. 

EPS2  -  Spread  criteria  for  multiple  roots.  (Input) 

ETA  If  the  root  X(I)  has  been  computed  and  it  is 

found  that  ABS (X(I)-X(J)) . LT. EPS2 ,  where  X(J)  is 
a  previously  computed,  then  the  computation  is 
restarted  with  a  guess  equal  to  X(I)+ETA. 

NSIG  -  Convergence  criterion.  (Input) 

A  root  is  accepted  if  two  successive 
approximations  to  a  given  root  agree  in  the 
first  NSIG  digits. 

N  -  The  number  of  roots  to  be  found.  (Input) 

X  -  Vector  of  length  N.  (Input/Output) 

On  input:  X  contains  the  initial  guesses  for  the 
roots . 

On  output:  X  contains  the  computed  roots. 

ITMAX  -  Iteration  indicator.  (Input/Output) 

On  input:  ITMAX  is  the  maximum  number  of 
iterations  to  be  taken  per  root. 

On  output:  ITMAX  is  the  number  of  iteration  used 
in  finding  the  last  root. 

IER  -  Error  parameter.  (Output) 

Warning  With  Fix: 

IER=33,  Indicates  that  for  one  root,  convergence 
was  not  obtained  within  ITMAX  iterations. 

That  root  is  set  to  111111. 

IER=34,  Indicates  that  for  one  root,  the 

derivative  of  the  function  at  that  root  was 
too  small.  That  root  is  set  to  222222. 
IER=35,  Indicates  that  the  error  conditions 
described  for  IERs33  and  IERS34  above 
occurred  more  than  once.  The  roots  fir  which 
the  error  occurred  are  set  to  111111  or 
222222,  depending  on  the  type  of  error. 
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PRECISION/HARDWARE  -  Single  and  Double/H32 

-  Single/H36,H48,H60 

REQD  IMSL  ROUTINES  -  UERTST ,  UQETIO 

NOTATION  -  Information  on  special  notation  and  conventions  is 

available  in  the  manual  introduction  or  through 
IMSL  routine  UHELP . 

REMARKS  1.  ZREAL2  assumes  that  there  exist  N  distinct  real  roots  for 

the  function  F  and  that  the  initial  guesses  supplied  by  the 
user  are  sufficiently  close  to  roots  to  obtain  convergence 
by  Newtons  method.  The  routine  is  designed  so  that 
convergence  to  any  single  root  cannot  be  obtained  from  two 
different  initial  guesses.  This  routine  is  intended 
primarily  for  the  refinement  of  N  known  rough 
approximations  of  the  roots  of  F. 

2.  Scaling  the  X  vector  in  the  function  F  may  be  required  if 
any  of  the  roots  are  known  to  be  less  than  one. 
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APPENDIX  C:  Data  Generation  Software 


SLAM  Code  for  Simulation  Model 


5  5  25  2000.  1 

100.0  0.0  1.0  2.78  2.78  25.0  25.0 

0.00  1.00  0.00  0.00  0.00  0.00  0.00 

0.00  0.00  1.00  0.00  0.00  0.00  0.00 

0.20  0.00  0.00  0.36  0.36  0.04  0.04 

0.00  0.00  1.00  0.00  0.00  0.00  0.00 

0.00  0.00  1.00  0.00  0.00  0.00  0.00 

0.00  0.00  1.00  0.00  0.00  0.00  0.00 

0.00  0.00  1.00  0.00  0.00  0.00  o.oo 

GEN, BAUER .DATA  GENERATION  MODEL , 06/05/86 , 1000 ,N ,N ,Y , N ,N; 
LIMITS, 7, 5, 200; 

STAT.l, RESPONSE  TIME; 

STAT , 2 .WAIT  STAT  2; 

STAT ,3 , WAIT  STAT  3; 

ST AT, 4, WAIT  STAT  4; 

STAT, 5, WAIT  STAT  5; 

STAT, 6, WAIT  STAT  6; 

STAT, 7, WAIT  STAT  7; 

TIMST.XXU)  .TERMINALS; 

TIMST ,XX(2) .CPU; 

TIMST ,XX(3) ,DISK1; 

TIMST, XX(4) .DISK2; 

TIMST. XX(5) .DISK3; 

TIMST, XX (6) .DISK4; 

TIMST, XX(7) .DISK5; 

INITIALIZE, 0. ,5000. ; 

MONTR, CLEAR, 2000; 

SEEDS, 34444866917(1) ; 

FIN; 
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o  o  o  o  o  —  o  o  o  o  o  o  o 


FORTRAN  Code  for  Simulation  Model 


PROGRAM  MAIN ( INPUT .OUTPUT ,TAPE5= INPUT ,TAPE6=00TPUT ,TAPE7 ,TAPE1 , 
&TAPE2 , TAPES , TAPE4 ) 

*t«#»*»*«*»#****#**##**#****»************** 

*  MAIN  PROGRAM  » 

«*#*#*#*##*####*»*##*****#**#**#*#********* 

PROGRAM  MAIN 
DIMENSION  NSET {5000) 

COMMON  QSET (5000) 

C0MM0N/SC0M1 /  ATRIB(IOO) ,DD(100) ,DDL(100) .DTNOW, I I ,MFA .MSTOP .NCLNR 
1 . NCRDR , NPRNT , NNRUN , NNSET , NT APE , SS ( 1 00 ) , SSL ( 100 ) , TNEXT , TNOW . XX ( 1 00 ) 
COMMON/UCOM1 /  DEPART (10) ,RMEAN(10) .P(IO.IO) ,SERVT( 10) ,EC0UNT(2) 
C0MM0N/UC0M2/  ISUBCAP , NUSSSN , NUMCUST , TCLEAR , NSTUDY 
EQUIVALENCE  (NSET ( 1 ) .QSET ( 1 ) ) 

NNSET=5000 
NCRDR=5 
NPRNT =6 
NTAPE=7 

READ  (NCRDR,*)  I SUBCAP, NUSSSN, NUMCUST, TCLEAR. NSTUDY 
READ  (NCRDR,*)  (RMEAN(I) , 1= 1 .NUSSSN+2) 

DO  10  1=1 .NUSSSN+2 

READ  (NCRDR,*)  (P ( I ,J) ,J=1 .NUSSSN+2) 

0  CONTINUE 

OPEN  (UNIT=10,FILE='DGM.0P' , STATUS = 'NEW’ ) 

CALL  SLAM 

CLOSE (10) 

STOP 

END 

i>mioHii»ii>«im*mi«*ii«»m*HH*H 

*  SUBROUTINE  EVENT  » 

•I***************************************** 

SUBROUTINE  EVENT ( I ) 

C0MM0N/SC0M1 /  ATRIB(IOO) ,DD( 100) ,DDL(100) , DTNOW , I I , MFA , MSTOP , NCLNR 
1 , NCRDR , NPRNT , NNRUN , NNSET , NT APE , SS ( 100) , SSL ( 100) , TNEXT , TNOW , XX ( 1 00 ) 
COMMON/ UC0M1 /  DEPART (10) , RMEAN(IO) ,P ( 10 , 10) ,SERVT(10) ,EC0UNT(2) 
C0MM0N/UC0M2 /  I SUBCAP , NUSSSN , NUMCUST , TCLEAR , NSTUDY 

ECOUNT ( 1 ) =ECOUNT ( 1 ) ♦ 1 

IF  (TNOW. GT. TCLEAR)  ECOUNT (2) =ECOUNT (2) +1 
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GOTO  (1,2)  ,  I 


1  CALL  ARSS 
RETURN 

2  CALL  ENDSS 
RETURN 

END 

«««•*«•**•*««»«««*•«««••««««»«« 

»  SUBROUTINE  INTLC  * 

##*#«#*«#«*§*«#***»#**#*****«***##»#**##*** 

SUBROUTINE  INTLC 

COMMON/SCOM1/  ATRIB(IOO)  ,DD<100)  .DDL(IOO)  ,DTNOW,II , MFA , MSTOP , NCLNR 
1 , NCRDR , NPRNT , NNRUN , NNSET , NT APE ,SS(100),SSL(100), TNEXT , TNOW , XX ( 1 00 ) 
C0MM0N/UC0M1/  DEPART (10) , RMEAN(IO) ,P( 10 . 10) ,SERVT(10) ,EC0UNT(2) 
COMMON /UC0M2/  ISUBCAP , NUSSSN , NUMCUST , TCLEAR , NSTUDY 
C0MM0N/GC0M5/  IISED(IO) , JJBEG , JJCLR , MMNIT ,MMON , NNAME (5) .NNCFI , 
&NNDAY , NNPT , NNPRJ ( 5 ) , NNRNS , NNSTR , NNYR , SSEED (10), LSEED (10) 
C0MM0N/UC0M3/  MULTIN0(7) 

INTEGER  ISEED (2000) 

IF  (NNRUN. EQ.l)  THEN 
DO  10  1=1,2000 

ISEED ( I ) = ( 1 . E+ 12) «DRAND ( 1 ) 

10  CONTINUE 

END  IF 

USED  (2)  =  ISEED  (NNRUN) 

X=DRAND (-2) 

DO  15  1=1,7 
MULTINO(I) =0 
15  CONTINUE 

DO  20  1=1,2 
ECOUNT ( I ) =0 . 

20  CONTINUE  " 

DO  25  I=l,NUSSSN+3 
DEPART (I) =0. 

25  CONTINUE 

DO  30  1=1 , NUSSSN+2 
SERVT (I ) =0 . 

30  CONTINUE 

DO  35  1=1, NUMCUST 

ETIME=EXPON(RMEAN( 1) ,2) 

ATRIB(1)=ETIME 
ATRIB (3) =1 
ATRIB(4) =1 
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ATRIB (5) =2 

CALL  SCHDL(1 ,ETI ME, ATRIB) 

35  CONTINUE 

DO  40  1=1 .NUSSSN+2 
XX(I)=0. 

40  CONTINUE 

WRITE (6, 100)  NNRUN.NSTUDY 

100  FORMAT (IX, ’SIMULATION  STUDY  IN  PROGRESS  :  RUN  ’,14,  *  OF 
&’  14 , '  RUNS’) 

RETURN 
END 

*##****##*#«**»####*»**«**#««»****«**«*»»»* 

*  SUBROUTINE  ENDSS  « 

a****#***##*#####*#*##**###****#*##*##*#**# 

SUBROUTINE  ENDSS 

COMMON/SCOM1/  ATRIB (100) ,DD(100) .DDL (100) , DTNOW , I I , MFA , MSTOP , NCLNR 
1 , NCRDR . NPRNT , NNRUN , NNSET , NTAPE , SS ( 1 00) , SSL ( 1 00) , TNEXT , TNOW , XX ( 1 00 ) 
C0MM0N/UC0M1/  DEPART (10) ,RMEAN(10) ,P(10,10) ,SERVT(10) ,EC0UNT(2) 
C0MM0N/UC0M2/  ISUBCAP , NUSSSN , NUMCUST , TCLEAR , NSTUDY 
C0MM0N/UC0M3/  MULTIN0(7) 

CALL  SCHDL( 1,0. .ATRIB) 

MYQ=ATRIB(4) 

IF  (NNQ(MYQ) .NE.O)  THEN 

CALL  RMOVEd  ,MYQ, ATRIB) 

WAIT=TN0W- ATRIB (2) 

CALL  COLCT( WAIT, MYQ) 

RM=RMEAN(MYQ) 

SERVICE=EXPON(RM, 2) 

ATRIB (4) = ATRIB (5) 

IAT=ATRIB (4) ♦ . 00001 
CALL  NEXTGUY ( I AT , I NEXT ) 

#*»##*«#»****»»***»**##»**#**»*#« 

»  COLLECT  STATISTICS  WHILE  » 

*  PARKED  AT  CPU  * 


IF  (IAT.EQ.3)  THEN 

MULT I NO ( I NEXT ) =  MULT I  NO  < I NEXT )  + 1 
END  IF 

ATRIB (5) "INEXT 

CALL  SCHDL( 2, SERVICE, ATRIB) 

IF  (TNOW. GT. TCLEAR)  THEN 

SERVT (MYQ) =SERVT (MYQ) ^SERVICE 
DEP ART ( MYQ )=  DEPART ( MYQ )  + 1 
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DEPART (NUSSSN+3) =DEPART (NUSSSN+3) +1 
END  IF 
ELSE 

XX(MYQ) =0 . 

END  IF 

IF  ( (MYQ.EQ.3) .AND. (NNQ(2) .GT.O) .AND. (ISUBCAP.NE.O) . AND. 
&(INEXT.EQ. 1) .AND. (NNQ(MYQ) .NE.O))  THEN 
CALL  RMOVE ( 1 , 2 , ATRIB) 

SERVICE=0 . 

ATRIB(4) =ATRIB (5) 

ATRIB (5) =3 

CALL  SCHDL(1 .SERVICE, ATRIB) 

END  IF 

RETURN 

END 

«  SUBROUTINE  ARSS  • 

*»###*####*#**#*##*******##**####***♦#»*»»* 

SUBROUTINE  ARSS 

C0MM0N/SC0M1/  ATRIB(IOO) ,DD(100) ,DDL(100) , DTNOW, I I , MFA , MSTOP , NCLNR 
1 , NCRDR . NPRNT , NNRUN , NNSET , NT APE , SS ( 100) , SSL ( 1 00) . TNEXT . TNOW, XX (100) 
C0MM0N/UC0M1/  DEPART (10) ,RMEAN(10) ,P(10,10) ,SERVT(10) ,EC0UNT(2) 
COMMON /UC0M2/  I SUBCAP , NUSSSN . NUMCUST , TCLEAR , NSTUDY 
C0MM0N/UC0M3/  MULTIN0(7) 

IAT=ATRIB(5) 

IF  (IAT.EQ.l)  THEN 
RESP=TNOW- ATRIB ( 1 ) 

CALL  COLCT(RESP.l) 

RM=RMEAN ( 1) 

SERVICE=EXPON (RM. 2) 

ATRIB ( 1 ) =TNOW+SERVICE 
ATRIB (4) =1 
ATRIB (5) =2 

CALL  SCHDL(1 .SERVICE .ATRIB) 

IF  (TNOW. GT. TCLEAR)  SERVT(IAT) =SERVT(IAT) + SERVICE 
GO  TO  101 
END  IF 

IF  (IAT.EQ.2)  THEN 

IF  (ISUBCAP.NE.O)  THEN 
NUMSUBsO 

DO  10  Is3 .NUSSSN+2 

NUMSUB=NUMSUB+NNQ (I) +XX(I) 

CONTINUE 

IF  (NUMSUB.LT. I SUBCAP)  THEN 
IF  (NNQ(2) . EQ.O)  THEN 
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WAIT=0 . 

CALL  COLCT (WAIT ,2) 

SERVICE=0 . 

ATRIB (4) =  2 
ATRIB (5) =3 

CALL  SCHDLd,  SERVICE,  ATRIB) 
GO  TO  101 
ELSE 

ATRIB (2) =TNOW 
CALL  FILEM(2, ATRIB) 

CALL  RMOVE (1,2, ATRIB) 
WAIT=TN0W-ATRIB(2) 

CALL  COLCT (WAIT, 2) 

ATRIB (4) =2 
ATRIB (5) =3 
SERVICED. 

CALL  SCHDL(1 .SERVICE, ATRIB) 
GO  TO  101 
END  IF 
ELSE 

ATRIB (2) =TNOW 
CALL  FILEM(2, ATRIB) 

RETURN 
END  IF 
END  IF 
END  IF 

100  IF  (XX(IAT) .GT.O. )  THEN 
ATRIB (2) =TNOW 
CALL  FILEMdAT .ATRIB) 

RETURN 

ELSE 

WAIT*0 . 

CALL  COLCT (WAIT, I AT) 

RM=RMEAN(IAT) 

ATRIB (4) =IAT 

CALL  NEXTGUY ( I AT , I NEXT ) 


«  COLLECT  STATISTICS  WHILE  * 
*  PARKED  AT  CPU  « 


IF  (IAT.EQ.3)  THEN 

MULT I NO ( INEXT) =MULTINO ( INEXT) + 1 
END  IF 

ATRIB (5) SINEXT 
SERVICE=EXPON (RM, 2) 

XX (I AT) = 1 

CALL  SCHDL( 2, SERVICE, ATRIB) 

IF  (TNOW . GT . TCLEAR)  SERVT(IAT) *SERVT(IAT) ♦SERVICE 
ENDIF 
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101  IF  (TNOW.GT. TCLEAR)  THEN 

DEPART ( I AT ) =  DEPART ( I AT ) + 1 
DEPART (NUSSSN+3) ^DEPART (NUSSSN+3) ♦ 1 
END  IF 

RETURN 
END 

*«*«•*«*««««**«««**«*»«*«••«»«*««***«»«««*« 

*  SUBROUTINE  NEXTGUY  * 

»##*###*******#******#*****»*»**#»*#**»»**« 

SUBROUTINE  NEXTGUY (I AT , I NEXT) 

COMMON/UCOM1/  DEPART (10) , RMEAN(IO) , P ( 10 , 10) ,SERVT(10) ,ECOUNT(2) 
COMMON/UCOM2/  ISUBCAP , NUSSSN , NUMCUST , TCLEAR , NSTUDY 

CUM=0 . 

U=UNFRM(0 . ,1. ,2) 

DO  10  INDEX* 1 , NUSSSN+2 
CUM=CUM+P (IAT , INDEX) 

IF  (U.LE.CUM)  THEN 
INEXT=INDEX 
GOTO  15 
ELSE 

CONTINUE 
END  IF 
CONTINUE 

RETURN 
END 

»  SUBROUTINE  OTPUT  » 


SUBROUTINE  OTPUT 

COMMON /SCOM1/  ATRIBdOO) ,DD( 100) ,DDL(100) .DTNOW.II , MFA , MSTOP , NCLNR 
1 , NCRDR , NPRNT , NNRUN , NNSET , NT APE , SS ( 1 00 ) , SSL ( 1 00 ) , TNEXT , TNOW , XX ( 1 00 ) 
COMMON/UCOM1 /  DEPART (10) , RMEAN(IO) ,P ( 10 . 10) ,SERVT(10) ,ECOUNT(2) 
COMMON/UCOM2/  ISUBCAP , NUSSSN , NUMCUST , TCLEAR , NSTUDY 
COMMON/UCOM3/  MULTIN0(7) 

WRITE ( 1 0 , » )  NNRUN 
WRITEdO,*)  (ECOUNT(I)  ,1*1,2) 

WRITEdO,*)  (CCAVG(I) , 1=1 .NUSSSN+2) 

WRITEdO,*)  (TTAVGd)  ,1=2, NUSSSN+2) 

WRITE(10,»)  (SERVT(I) , 1= 1 .NUSSSN+2) 

WRITE(10,»)  (DEPART ( I ) ,1=1, NUSSSN+3) 
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ISUM=0 
DO  10  1=1,7 

I SUM= I SUM* MULT I NO ( I ) 

10  CONTINUE 

WRITE (10,*)  (MULTINO(I)  ,1  =  1,7)  ,ISUM 

RETURN 

END 
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FORTRAN  Code  for  Data  Post-Processing 


PROGRAM  MKDATA 
C 

C  tHHmHHHItHHtUKHMHfMMHMHMMHtHMIHUittHHtH 
C  «  THIS  PROGRAM  TAKES  THE  VARIABLE  OUTPUTS  FROM  THE  DATA  GENERATION  « 
C  *  MODEL  AND  COVERTS  THEM  TO  WORK/ROUTING  VARIABLES  AND  RESPONSES.  • 
C  *  THESE  VARIABLES  POSSESS  THE  DESIREABLE  QUALITIES  OF  A  MEAN  OF  * 
C  *  ZERO  AND  VARIANCE  OF  ONE.  THIS  ‘CONVERTED*  DATA  IS  THEN  USED  AS  • 
C  *  INPUT  TO  THE  VARIABLE  SUBSET  SELECTION  PROGRAM.  • 

C  #«««****##«#«****###«*##**«**##***»#*«#*#**«»**»»»#**»**»»»*#»»»»**» 
C 

PARAMETER  (NUMREPS8 1000) 

C 

REAL  R( 13) ,W(7) ,E(8) 

REAL  WRKV (7) 

REAL  RM(7) ,PI(7) ,YM(13) 

REAL  VEC(4) ,RMULT(7) ,PI2(7) 

REAL  WK (7)  ,  WK2  (7) 

INTEGER  MULT (8) ,DP 
C 

C  *»******#**#*«#*««*#*##***«*»#i»#*##*«*#»*»«*******#***********»«»»»* 

C  *  DATA  STATEMENTS  » 

C  «  « 

C  »  VEC  8  VECTOR  OF  STEADY  STATE  RESPONSE  MEANS  * 

C  «  RM  8  VECTOR  OF  MEAN  SERVICE  TIMES,  BY  STATION  * 

C  »  PI  =  VECTOR  OF  STEADY  STATE  TRANSITION  PROBABILITIES  * 

C  «  (DERIVED  ANALYTICALLY)  * 

C  »  PI2  =  VECTOR  OF  ACTUAL  BRANCHING  PROBABILITIES  FROM  CPU  • 

C  »  (AS  SUPPLIED  TO  THE  SLAM  PROGRAM)  « 

C 

DATA  VEC  /30.72,  1.047,  1.458,  13.09/ 

DATA  RM  / 100.0,  0.0,  1.0,  2.78,  2.78,  25.0,  25.0/ 

DATA  PI  /0.00,  0.09,  0.45,  0.16,  0.16,  0.02,  0.02/ 

DATA  PI2  /0 . 2 ,  0.0,  0.0,  0.36,  0.36,  0.04,  0.04/ 


«  OPEN  INPUT  AND  OUTPUT  FILES  * 


OPEN  (UNIT8 10 ,  FILE8 ' DGM. OP ’ ,  STATUS8 ’OLD’ ) 

OPEN  (UNIT=20 ,  FILE='VSSP1.IN' ,  STATUS8 ' NEW' ) 

OPEN  (UNIT*30 ,  FILE8’ SUMMARY’ ,  STATUS8 ’NEW* ) 

C 

C  «  READ  DATA  FROM  FILES  AND  CONVERT  TO  WORK  AND  ROUTING  VARIABLES 
C  * 

C  •  DATA: 

C  *  R1  8  RUN  • 

C  *  R2  8  EVENT  COUNT 


C  •  R3  *  EVENT  COUNT  FROM  TCLEAR  » 

C  «  R  =  VECTOR  OF  RESPONSES  * 

C  »  W  =  VECTOR  OF  SUMS  OF  SERVICE  TIMES,  BY  STATION  * 

C  *  E  =  VECTOR  OF  TOTAL  DEPARTURES,  BY  STATION  * 

C  *  MULT { J)  =  VECTOR  OF  TOTAL  DEPARTURES  FROM  CPU  TO  STATION  J  * 

C  *  « 

C  *  OUTPUT:  « 

C  «  WRKV(J)  =  WORK  VARIABLE  J  * 

C  «  RMULT ( J )  =  ROUTING  VARIABLE  J  » 

C  ***#******«##*«***#***#*«**##*»#«*#***»#**«»#»»»t»#»»#»*#**#*»»***»* 

c 

DO  10  1=1, NUMREPS 
C 

READ ( 10 , »)  R1 

READ (10,*)  R2,R3 

READ ( 1 0 , « )  (R(II) ,11=1 ,7) 

READ ( 1 0 , « )  (R(II) ,11=8,13) 

READ (10 , «)  (W(II) ,11=1,7) 

READ(10,*)  (E(II) ,11=1 ,8) 

READ ( 10 , «)  (MULT(II) ,11=1 ,8) 

C 

DO  15  J= 1 ,7 
C 

IF  (RM(J) .NE.O. )THEN 

WRKV (J)=(W(J)-E(J) *RM( J) ) « (SORT (E ( J) ) 

&  / (PI ( J) *E (8) «RM(J) ) ) 

END  IF 
C 

IF  ((J.EQ.2) .OR. (J.EQ.3))  THEN 
RMULT (J)=0. 

ELSE 

RMULT  (J)  =  (MULT(J)  -MULT  (8)  »PI2  (J) ) 

&  /SORT (PI2(J) « (1 PI2(J) ) »MULT (8) ) 

END  IF 
C 

15  CONTINUE 

C 

DO  20  J= 1 , 13 

YM( J) =YM(J) +R(J) 

20  CONTINUE 

C 

DO  25  J=1 ,7 

WK(J)=WK(J)+WRKV(J) 

WK2 (J) =WK2 (J) +WRKV (J) *«2 
25  CONTINUE 

C 

C  WRITE(20 ,»)  RMULT(l) ,RMULT(4) ,RMULT(5) ,RMULT(ti) ,RMULT(7) 

WRITE (20 , «)  WRKV ( 1 ) ,WRKV(4) ,WRKV(5) ,WRKV(6) ,WRKV(7) , 

&  R(l) ,R(9) , R( 10) ,R ( 1 1) ,R( 12) ,R(13) 

C 

10  CONTINUE 
C 
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C  «  CALCULATE  THE  MEANS  OF  THE  RESPONSES  (YM) ,  AND  THE  MEANS  (WK)  * 

C  *  AND  VARIANCES  (WK2)  OF  THE  WORK  VARIABLES.  THEN  PRINT  SUMMARY  « 
C  *  INFORMATION  TO  A  SUMMARY  FILE.  * 

C  #*#*»»******»»#**»»»«»»*#»#**#*#»#»*#*«*«»*****»*»**»»**»**»******«« 

c 

DO  30  J=1 , 13 

YM(J)=YM(J)/ (FLOAT (NUMREPS) ) 

30  CONTINUE 

C 

DO  35  J=1 ,7 

WK ( J) =WX ( J) / (FLOAT (NUMREPS) ) 

WK2 ( J) =WK2 ( J) / (FLOAT (NUMREPS) ) 

35  CONTINUE 

C 

DO  40  J=1 ,7 

WK2 (J) =WK2 (J) -WK(J) »*2 
40  CONTINUE 

C 

WRITE (30, 555) 

WRITE (30, 556)  NUMREPS 
WRITE (30 , 557) 

WRITE (30, «)  VEC 
WRITE (30, 558) 

WRITE(30,«)  YM 
WRITE (30, 559) 

WRITE (30 , »)  WK 
WRITE (30, 560) 

WRITE (30 , *)  WK2 
C 

555  FORMAT (IX, ’SUMMARY  FILE  OF  DGM.OP  DATA  POST  PROCESSING’/) 

556  FORMAT ( IX, ’VSSP1. IN  HAS  A  TOTAL  OF  ’,15,’  REPLICATIONS’/) 

557  FORMAT (IX, ’BELOW  ARE  THE  POPULATION  MEANS  OF  THE  RESPONSES’/) 

558  FORMAT (IX, ’BELOW  ARE  THE  SAMPLE  MEANS  OF  THE  RESPONSES’/) 

559  FORMAT (IX, ’BELOW  ARE  THE  MEANS  OF  THE  WORK  VARIABLES’/) 

560  FORMAT (IX, 'BELOW  ARE  THE  VARIANCES  OF  THE  WORK  VARIABLES’/) 

C 

CLOSE (10) 

CLOSE (20) 

CLOSE (30) 

C 

STOP 

END 
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